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Preface 


In  early  1992,  Dr.  W.  E.  Wiesel,  a  professor  at  the  Air  Force  Institute  of  Tech¬ 
nology,  wrote  a  paper  postulating  the  existence  of  an  extended  Lyapunov  exponent.  He 
questioned  the  fact  that  the  classical  Lyapunov  exponent  contained  only  a  real  part,  while 
the  analogous  exponents  for  the  constant  coefficient  and  periodic  coefficient  cases  contain 
both  a  real  and  imaginary  part.  For  this  thesis,  I  developed  the  computer  code  to  calculate 
the  extended  Lyapunov  exponent,  as  defined  by  Wiesel,  for  .second  order  systems.  I  then 
applied  this  to  some  classic  trial  cases.  I  obtained  excellent  results  for  both  the  constant 
coefficient  and  periodic  coefficient  cases.  Thus,  the  technique  presented  here  validates 
the  existence  of  the  imaginary  part  of  the  Lyapunov  exponent.  Since  the  concept  of  an 
extended  Lyapunov  exponent  unifies  linear  system  theory,  this  work  should  definitely  be 
continued.  In  particular,  the  technique  needs  to  be  applied  to  a  wider  variety  of  systems. 

In  writing  the  code  and  analyzing  the  data,  I  received  a  great  deal  of  help  from 
others.  I  am  deeply  indebted  to  my  faculty  advisor,  Dr.  W.  E.  Wiesel,  for  his  untiring 
patience  and  assistance.  Additionally,  I  wish  to  thank  Captain  Chris  Hall  for  the  numerous 
impromptu  discussions  about  the  general  subject.  A  word  of  thanks  also  must  be  sent 
to  Major  Robinson,  Dr.  Phil  Beran  and  Mr.  Ron  Gordon  of  the  ENY  Computational 
Dynamics  and  Design  Laboratory  for  the  use  of  the  Sun  workstations  and  the  tailoring  of 
the  software  to  meet  my  specific  needs. 

I  must  also  recognize  the  endless  support  that  my  parents  have  provided  throughout 
my  life.  Thank  you  Mom  for  instilling  in  me  a  love  of  learning.  I  wish  you  could  be  here  to 
see  this,  though  I  know  you  are  here  in  spirit.  Thank  you  Dad  for  teaching  me  the  value 
of  stick-to-itiveness.  I’m  so  happy  that  you  can  share  this  with  me. 

Additionally,  I  wish  to  thank  my  wonderful  family.  To  Fred  and  Jesse:  Thank  you, 
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Abstract 

It  is  shown  that  the  concept  of  the  Lyapunov  exponent  can  be  extended  to 
include  an  imaginary  part.  A  numerical  technique  used  to  calculate  these  extended 
Lyapunov  exponents  for  second  order  systems  is  presented.  There  are  two  require¬ 
ments  to  make  this  extension  possible;  the  definition  of  a  coordinate  frame  on  the 
tangent  space  of  the  differential  equation,  and  an  extension  of  the  classical  limit, 
called  the  limit  of  the  mean.  An  application  of  the  technique  to  the  van  der  Pol 
equation  for  the  constant  coefficient  and  periodic  coefficient  cases  is  given.  The 
extended  Lyapunov  exponents  found  using  this  technique  totally  agree  with  the 
eigenvalues  for  the  constant  coefficient  case.  In  the  periodic  coefficient  case,  not 
only  do  the  extended  Lyapunov  exponents  agree  with  the  Poincare  exponents  calcu¬ 
lated  using  standard  Floquet  theory,  they  confirm  that  the  imaginary  parts  of  the 
Poincare  exponents  are  equal  to  the  quotient  where  t  is  the  period  of  the 

trajectorj’.  This  imaginary  part  is  uncertain  when  using  standard  Floquet  theory. 
Additionally,  fast  Fourier  transform  (FFT)  techniques  are  used  to  validate  the  exis¬ 
tence  of  the  extended  Lyapunov  exponent  and  the  values  obtained  for  its  imaginary 
part.  These  techniques  show  that  the  power  spectrum  of  relative  motion  is  discrete 
for  the  trial  cases  presented,  with  the  fundamental  frequency  almost  exactly  equal 
to  ihe  calculated  imaginary  part  of  the  extended  Lyapunov  exponent.  Coupled  with 
the  successful  comparison  of  characteristic  exponents  for  the  constant  coefficient 
and  periodic  coefficient  cases,  this  power  spectrum  serves  to  decisively  validate  the 
existence  of  the  extended  Lj  ipunov  exponent. 
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EXTENDED  LYAPUNOV  EXPONENTS 
FOR  SECOND  ORDER  SYSTEMS 


I.  Introduction 

Why  investigate  an  extended  Lyapunov  exponent?  In  virtually  every  area  of 
the  sciences,  problems  of  interest  involve  system  elemental  rates  of  change  dependent 
upon  the  interaction  of  the  system  basic  elements.  This  interaction  is  universally 
expressed  as  a  system  of  first  order  differential  equations 

X(t)  =  F(X,t)  (1.1) 

This  system  of  equations  could  model  anything  from  a  mechanical  system  to  n  chem¬ 
ical  reaction  or  an  economic  event.  Once  a  particular  solution,  Xp(t),  is  obtained 
for  1.1,  the  analyst  must  then  ask  the  following  questions: 

1)  How  stable  is  this  solution  ?  (i.e.  Will  some  slight  variation  cause  the 
system  to  diverge  from  the  particular  solution  ?) 

2)  VVdiat  is  the  system’s  fundamental  frec(uency:  at  what  frequency  does  it 
possess  energy  ? 

Obviously,  the  question  of  stability  is  paramount  in  predicting  system  behavior, 
but  the  need  to  understand  a  system's  driving  frequencies  is  equally  important. 
Improper  calculation  of  a  system’s  frequency  components  has  been  the  ultimate 
cause  of  numerous  disasters,  such  as  the  collapse  of  a  suspension  bridge  (8).  When  a 
suspension  bridge  was  excited  at  just  the  right  frecjuency,  the  extra  energy  injected 
into  the  system,  in  addition  to  the  energy  it  already  possesssed  at  that  frequency, 
caused  a  total  collapse  of  the  structure. 
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The  questions  of  stability  and  frequency  content  are  both  answered  by  a  sys¬ 
tem’s  characteristic  exponents,  where  the  real  parts  of  the  exponents  are  the  key  to 
the  stability  of  the  solution,  and  the  imaginary  parts  determine  the  system's  funda¬ 
mental  frequency.  The  eigenvalues  are  the  characteristic  exponents  for  a  constant 
coefficient  system,  and  the  Poincare  exponents  are  the  characteristic  exponents  for  a 
periodic  coefficient  system.  Similarly,  the  Lyapunov  exponents  are  the  characteristic 
exponents  for  a  general  time-varying  sytsem;  however,  classic  Lyapunov  exponents 
are  defined  to  be  real,  not  complex  numbers.  Thus,  they  yield  no  freqi  ency  infor¬ 
mation.  riie  extendtd  Lyapunov  exponents  disciussed  in  this  thesis  contain  both  real 
and  iiiinfjinary  parts.  The  following  section  outlines  the  classical  analysis  methods 
I  hat  were  used  as  the  basis  for  the  concept  of  extending  t  he  Lyapunov  exponent  to 
include  an  imaginary  part. 

/./  Classical  Analysis 

To  analyze  the  particular  .solution,  one  first  writes  a  nearby  trajectory  as 

X(t)-Xp(t)-f6X(t)  (1.2) 

where  (“iXft)  is  the  first  order  difference  l)etwcen  adjacent  trajectories  and  th<'  par¬ 
ticular,  or  known,  trajectory.  Substituting  equation  1.2  into  l.I  yields 

X(t )  =  Xp(t)  -h  6X(i)  =  F(Xp  -f  <sx.  t)  ( I.;i) 

Kx))anding  F  in  a  f’aylor  series  about  Xp,  and  assuming  i^X(t)  is  small,  produces 

X„(t)  +  ,^X(.)^F(Xp.t)+  ^  <SX  (1.1) 

f/X  V- 
Ap 


1-2 


Since  Xp(t)  =  F(Xp,t),  we  may  cancel  it  from  both  sides  of  1.4  to  obtain  the 
equation  of  variation 

^X(t)=^  SX{t)  (1.5) 

Xp 

rhe  partial  derivative  of  the  vector  F  with  respect  to  X  is  a  square  matrix  generally 
referred  to  as  a  Jacobian  matrix,  that  will  hereafter  be  referred  to  as  A(t).  In  essence, 
A(t)  is  the  matrix  of  partial  derivatives  of  the  differential  equations  with  respect  to 
the  variables  of  the  problem,  the  elements  of  X.  After  evaluation  on  the  particular, 
or  nominal,  trajectory  Xp,  A  is  a  function  of  time  alone. 

Since  equation  1.5  has  had  quadratic  terms  truncated  (due  to  small  <5X),  the 
character  of  its  solutions  governs  the  linear'  stability  of  the  trajectory  Xp(t).  Also, 
sincp  equation  1.4  is  a  linear  aystem,  its  solution  can  be  written  as 

<5X(t)  =  $(t,to)6X(to)  (1.6) 

where  the  fundamental  matrix  $  obeys  the  following 


^(t,to) 

=  A(t)$(t,to) 

(1.7) 

toi  to) 

=  I  (the  identity  matrix) 

(1.8) 

When  Xp  is  an  equilibrium  point  of  1.1,  then  A,  the  associated  linear  sys¬ 
tem  of  I.l.  contains  constant  terms.  For  linear,  constant  coefficient  systems,  the 
fundamental  matrix  can  be  written  as 

$(t.t„)  =  Eexp(J(t-to))E-‘  (1.9) 

where  the  eigenvectors  of  A  form  the  E  matrix,  and  the  J  matrix,  a  matrix  of  con¬ 
stants.  is  the  Jordan  normal  form  of  A  with  the  eigenvalues  of  A  along  its  diagonal. 
I'he  real  part  of  the  eigc  values  determine  stability;  n  gative  real  |)arts  imply  stabil- 


ity,  positive  real  parts  indicate  instability,  and  real  parts  equal  to  zero  imply  marginal 
stability  (8).  The  imaginary  parts  of  the  eigenvalues  determine  the  fundamental  or 
oscillation  frequencies. 

When  Xp  is  a  periodic  orbit,  A  is  a  periodic  matrix,  and  Floquet  theory  is 
applicable  (2).  In  this  instance,  the  fundamental  matrix  is  written  as 

$(t,to)  =  E(t)exp(J(t  -  to))E“'(to)  (1-10) 

Now,  however,  the  eigenvector  matrix  is  a  function  of  time,  and  is  periodic.  The 
diagonal  elements  of  J  are  now  called  Poincare  exponents.  Due  to  its  periodicity, 
E(t)  is  bounded.  Thus,  analogous  to  the  eigenvalues  in  the  constant  coefficient  case, 
all  stability  information  is  found  in  the  Poincare  exponents. 

When  Xp  is  a  general  trajectory,  A(l)  is  time  varying.  For  this  general  case, 
Lyapunov  (6)  showed  that  there  are  N  real  numbers,  Lyapunov  exponents,  which 
determine  the  local  exponential  rates  of  growth  away  from,  or  convergence  to  the 
reference  trajectory.  Accepted  algorithms  for  calculating  Lyapunov  exponents  have 
been  reported  by  Benettin,  et  al  (3)  and  by  Shimada  and  Nagashimi  (9).  Extending 
this  general  case  exponent  to  include  an  imaginary  part  is  the  heart  of  this  thesis. 

1.2  Problem  Statement 

As  stated  in  the  previous  section.  Lyapunov  showed  the  existence  of  N  real. 
not  complex,  numbers  characterizing  the  growth  away  from,  or  convergence  to  a 
particular  solution.  In  his  paper  “Extended  Lyapunov  Exponents'’(  1 1 ),  Wiesel  posed 
the  following  question:  Given  the  fact  that  in  the  constant  coefficient  and  periodic 
coefficient  cases,  the  characteristic  exponents  are  complex  numbers,  why  are  the 
Lyapunov  exponents  limited  to  being  real  ? 


Thus,  in  light  of  1.9  and  1.10,  Wiesel  proposed  a  fundamental  matrix  decom¬ 
position  for  the  general  case  to  be 

^(t,to)  =  E(t)exp(J(t  -  to))E“^(t)  (1-11) 

where  the  eigenvector  matrix  E(t)  would  take  on  the  character  of  A(t):  constant, 
periodic,  multiply  periodic,  or  aperiodic,  and  the  diagonal  elements  of  the  constant 
Jordan  matrix  J  would  govern  stability. 

The  main  objectives  of  this  study  are  to  calculate,  and  validate  the  existence 
of  extended  Lyapunov  exponents  (i.e.  ones  that  include  an  imaginary  part),  for  some 
trial  cases.  The  existence  of  an  extended  Lyapunov  exponent  would  unify  the  theory 
of  linear  systems  for  all  cases.  Computer  programs  need  to  be  developed  to  apply 
the  technique  to  some  trial  cases.  During  this  study,  the  trial  cases  will  be  limited  to 
second-order  systems  .  After  development  of  the  computer  programs,  the  next  step 
will  be  to  apply  them  to  a  second  order  system  and  analyze  some  classic  behaviors: 
the  constant  coefficient  and  periodic  coefficient  cases.  Comparison  of  the  calculated 
extended  Lyapunov  exponents  to  the  eigenvalues  of  the  A  matrix  for  linearization 
about  an  equilibrium  point  will  validate  the  constant  coefficient  case.  Linearizing 
about  the  limit  cycle,  provided  a  limit  cycle  exists,  will  serve  to  validate  the  periodic 
coefficient  case.  Naturally,  for  the  periodic  coefficient  case,  the  comparison  will  be 
against  the  Poincare  exponents  which  are  calculated  using  Floquet  theory. 

1.3  Overview  of  Thesis  Content 

The  following  chapter  discusses  the  derivation  of  the  extended  Lyapunov  ex¬ 
ponent.  It  includes  both  an  heuristic  derivation  and  a  derivation  based  on  Lya¬ 
punov's  original  definition.  Additionally,  it  covers  the  development  of  the  limit  of 
the  mean  and  the  coordinate  frame  definition.  Chapter  III  contains  a  brief  overview 


1-.5 


of  the  numerical  algorithm  development  and  the  associated  FORTRAN-coded  pro¬ 
grams/subroutines. 

In  Chapter  IV,  the  reader  will  find  an  application  of  the  technique  to  a  second 
order  system,  van  der  Pol’s  equation,  with  linearization  about  the  equilibrium  point 
and  the  limit  cycle.  Also,  Chapter  IV  contains  the  FFT  analysis  used  to  validate 
the  imaginary  parts  of  the  extended  Lyapunov  exponents.  The  last  chapter  contains 
the  author’s  conclusions  and  recommendations. 


1-6 


II.  Derivation  of  the  Extended  Lyapunov  Exponent 


Wiesel’s  (11)  arguments  for  extending  the  Lyapunov  exponent  will  be  presented 
in  this  chapter.  The  first  section  will  contain  two  derivations  of  the  extended  Lya¬ 
punov  exponents,  the  first  heuristic  and  the  second  generalizing  Lyapunov's  original 
definition.  It  will  be  shown  in  the  second  section  that  the  definition  of  a  coordinate 
frame  is  essential  to  the  calculation  of  an  imaginary  part  of  an  exponent,  even  in  the 
constant  coefficient  case.  An  extended  definition  of  a  limit,  the  limit  of  the  mean, 
will  be  developed  in  the  third  section.  Finally,  the  last  section  imparts  VViesel’s 
original  conjecture. 

2.1  The  Fundamental  Matrix  and  the  Extended  Lyapunov  Exponent 
A  basic  property  of  the  fundamental  matrix  $  is 

$(t,to)  =  A(t)#(t,to)  (2.1) 

where  t  is  the  current  time,  and  to  is  the  initial  time.  Equation  2.1  is  a  general, 
time  dependent  linear  system  with  the  initial  condition,  $(to,to)  =  I  (the  identity 
matrix).  As  will  be  shown  in  the  forthcoming  section,  this  initial  condition  is  crucial. 

Using  standard  matrix  eigenvector  decomposition 

$(t,to)  =  E(t)A(t)E-’(t)  (2.2) 

where  E(t)  and  A(t)  are  the  fundamental  matrix  eigenvector  and  eigenvalue  matrices, 
respectively.  It  is  assumed  that  the  eigenvalues  A,  on  the  diagonal  of  the  diagonal 
matrix  A  are  all  distinct,  so  that  A  is  not  a  Jordan  form.  Since  4^(to,to)  =  I,  A(to) 
must  equal  the  identity  matrix,  and  E(to)  should  be  .some  invertible  matrix.  If  the 
eigenvectors  are  normalized,  then  any  long  term  stability  information  resides  within 
A. 
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To  further  continue  the  analogy  with  1.11,  a  diagonal  matrix  17  is  defined  by 


—  to)  =  In  Ai(t)  (2-3) 

where  the  initial  condition  A(to)  =  I  is  satisfied.  Equation  2.2  can  now  be  rewritten 
as 

$(t,  to)  =  E(t)exp(17(t)(t  -  to))E-'(t)  (2.4) 

Comparison  of  equation  2.4  to  1.11  shows  that  the  exponential  rates  of  change 
of  the  solution  vectors  over  the  interval  (to,t)  are  given  by  u;,(t),  where 

u^.(t)  =  ^^^!nAi(t)  (2.5) 

Then,  should  the  limit  exist  as  time  approaches  infinity,  the  long-term  system  re¬ 
sponse  is  defined  by 

u;,(cxd)  =  limu;,(t)  (2.6) 

t—^00 

where  the  real  parts  of  w,(t)  are  the  average  exponential  rates  of  increase  of  the  solu¬ 
tion,  and  the  imaginary  parts  are  the  average  oscillation  frequencies  of  the  solution 
with  respect  to  the  eigenvector  basis.  (The  eigenvectors  are  independent,  and  span 
the  .solution  vector  space.) 

Thus,  the  real  parts  are  the  classic  Lyapunov  exponents  of  the  system.  The 
imaginary  parts  must  exist  in  comparison  to  the  constant  coefficient  and  periodic 
coefficient  cases,  though  no  mention  of  them  was  found  in  any  literature  searches 
conducted.  Therefore,  Wiesel  dubbed  u;,(oo)  as  the  extended  Lyapunov  exponents. 

The  following  derivation  will  use  a  different  tack  to  arrive  at  equation  2.6: 
generalizing  Lyapunov’s  original  definition.  First,  taking  equation  2.4,  and  post 
multiplying  the  left  and  right  sides  by  E(t)  yields 

$(t,to)E(t)  =  E{t)exp[J7(t)(t  -  to)]  (2.7) 
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or, 

ei(t)exp[aJi(t)(t  -  to)]  =  4>(t,  to)ei(t)  (2.8) 

where  ei(t)  are  columns  of  E(t).  The  revealing  form  of  equation  2.8,  upon  comparison 
to  1.6,  shows  that  each  column  vector  of  equation  2.8  is  a  solution  to 

lii(t)  =  A(t)ui(t)  (2.9) 

for  the  following  initial  conditions  and  final  values: 

Ui(to)  =  ei(t)  (2.10) 


Ui(t)  =  ei(t)exp(a;i(t)(t  -  to))  (2.11) 

It  is  assumed  that  Ui(to)  is  normalized. 

Next,  Lyapunov  (6)  showed  that  for  iq(t)  =  A(t)ui(t),  N  real  Lyapunov  expo¬ 
nents  exist,  which  obey 


fJ-i 


lim  sup  ki(t)exp(u;i(t)(t-to))| 

(t-to) 

lini  -I-  |exp(u;i(t)(t-to))| 

t-oc  f^[(t-to)  (t-to) 


(2.15) 

(2.16) 


2- .3 


Since  the  eigenvectors  were  assumed  to  be  normalized,  the  following  statement  is 
true: 


In  |ei(t)|  =  0 


(2.17) 


Thus, 


//,  =  lim  sup 


In  |exp(u;i(t)(t  -  to))| 


t-*oo  ‘  (t  -  to) 


Now,  let  u;,(t)  =  0  +  1/3,  where  i  =  \/— T. 


Hi  =  lim  sup 

t—^OO 


In  |exp[(Q  +  i/j)(t  -  to)]  I 
(t  -  to) 

ln[[|exp(o(t  -  to))|][|exp(i^(t  -  to 


lim  sup 

*-*'>=  (t-to) 


=  lim  sup 


In  |exp(o(t  -  to))j  ^  In  |exp(i/3(t  -  to))| 


(t  -  to) 


(t-to) 


(2.18) 


(2.19) 

(2.20) 
(2.21) 


Making  use  of  Euler’s  formula  (4):  e‘®  =  cos(0)  +  isin(^). 


Hi  =  lim 

*OCi 


In  |exp(Q(t  -  to))|  ,  In  |cos(/3(t  -  to))  +  isin(/3(t  -  to) 
sup - - — — - +  sup - — 


(t  -  to) 


(t  -  to) 


(2.22) 


For  suprema  conditions,  the  second  In  term  vanishes.  Thus, 


Hi  =  lim  sup 


In  |exp(Q(t  —  to) 


‘  (t  -  to) 


lim  sup(q) 


(2.23) 


=  lim  sup /?e  a;,(t) 

*-oo 


Since  the  eigenvectors  are  independent  and  span  the  solution  space  (i.e.  form 
a  basis),  the  N  trial  trajectories  span  the  space  of  initial  conditions,  and  any  other 
trajectory  with  |u(to)|  =  1  can  be  written  as  a  linear  combination  of  Ui(to).  This 
implies  that  the  classical  Lyapunov  exponents  are  (provided  the  limit  exists) 

//,  =  lim  .sup  77-^-— /fe  In  A,(t)  (2.24) 

(t  —  toj 
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where  A,  are  the  eigenvalues  of 

It  is  highly  unlikely  that  Ui(t)  and  Ui(to)  will  be  parallel  in  the  classic  case,  thus 
the  norm  is  required.  Since  the  product  of  a  matrix  and  a  vector  equals  a  vector, 
and  the  norm  of  a  vector  is  equal  to  a  real  number,  it  is  obvious  from  equation 
2.12  that  it  is  this  norm  that  limits  the  classical  Lyapunov  exponents  to  being  real, 
rather  than  complex,  numbers.  To  show  that  this  is  true,  the  reader  should  note  the 
following  definitions  for  the  norm  of  a  vector. 

First,  if  the  vector  y  is  a  vector  of  real  numbers,  its  norm  is  defined  (10) 

|y|  =  \lyl  +  y2  +  +  (2-25) 


where  t/,  are  the  elements  of  y,  and  the  norm  of  y  is  a  real  number.  Next,  if  y 
contains  complex  numbers,  the  following  definition  for  its  norm  applies  (10) 


|y|  =  yjyxyi  +*/2y2  +  •••  + 


,Vny„ 


(2.26) 


where  IJ-  is  the  complex  conjugate  of  y,.  Thus,  since  the  product  of  a  complex  number 
and  its  conjugate  is  always  a  positive  real  number,  the  norm  of  the  complex  vector 
y  is  real. 

Now,  since  this  derivation  has  used  the  eigenvectors  as  trial  trajectories,  the 
norm  becomes  unnecessary.  This  fact  is  highlighted  by  equation  2.8.  Which  shows 
that  when  propogating  ui  from  to  to  t,  Ui(to)  =  ei(t)  is  simply  multiplied  by  a  scalar, 
'fherefore,  Ui(t)  is  parallel  to  Ui(to),  which  eliminates  the  need  for  a  norm. 

Elimination  of  the  norm  from  2.12  yields 


Hi 


li^j  gy  In  [exp(u;i(t)(t  -  tp))] 
(t-to) 

=  lim  u:i{t) 

*OC 


(2.27) 


=  u.’,(oo) 


25 


Equation  2.27  is  identical  to  the  heuristically  derived  2.6. 


2.2  The  Need  for  Coordinates 

As  the  previous  section  illustrated,  the  definition  of  the  classical  Lyapunov 
exponents  assumes  that  the  tangent  space  to  the  differential  equation  X  =  F(X,  t) 
possesses  a  norm,  but  not  a  coordinate  frame.  This  section  will  use  two  simple 
examples  to  illustrate  the  crucial  role  that  the  choice  of  a  coordinate  frame  plays  in 
the  extended  Lyapunov  exponent. 

h'irst,  consider  the  simple  autonomous  linear  system,  where 


A  = 


Two  solutions  to  $  =  are 

^i(t) 

^2(t) 


(2.28) 


(2.29) 


In  either  solution,  one  column  vector  grows  with  logarithmic  rate  of  +1,  and 
the  other  decreases  with  logarithmic  rate  —1.  Clearly  then,  the  Lyapunov  exponents, 
which  are  also  the  eigenvalues  of  A,  are  ±1  in  either  case. 

The  next  step  is  to  calculate  the  extended  Lyapunov  exponents  for  the  two 
solutions.  First,  the  calculation  of  the  eigenvalues  of  ^  must  be  accomplished.  Fo 
do  this,  one  must  solve  the  characteristic  equation  which  is  defined  by 


det  {$  -  AI}  = 


<^]l  —  A  ^12 
<^21  *p22  ~  A 


=  0 


(2.30) 
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For  $1,  this  determinant  yields 


+  A(-e-‘  -e‘)  +  1  =  0 


(2.31) 


=>  A, 


'(e  ‘  +  e*)  ±  \J{e  ^ 

2 

(e“‘  +  e‘)  ±  (e“‘  —  e‘) 

2 


Thus,  Ai  =  e  ‘  and  A2  =  e‘. 

Solving  for  the  extended  Lyapunov  exponents  yields 


oj\  =  lim  -Ine*^  =  lim  -  =  I 

t-.oo  t  t 

1 ,  -t 

ijj2  =  lim  -Ine  =  lim  —  =  —  1 

t-.co  t  t— 00  t 

(i.e.  the  expected  values). 

Next,  the  characteristic  equation  for  ^2  is 


(2..32) 


(2.33) 

(2.34) 


A^  +  1=0  (2.35) 

=»  A,  =  ±i  (2.36) 

Another  Euler  identity  states  that  the  polar  form  of  a  complex  number  o  +  i/?  is 
expressed  as  re'^  where  r  =  sqrt(Q^  +  j5^)  and  9  =  arctan  |f|- 


=>•  Ai  =  +i  =  e"^  and  A2  =  — i  =  e 


(2.37) 
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Using  the  polar  form  of  A,  to  calculate  the  extended  Lyapunov  equations  yields  the 
following: 


lim  -Ine^'"^  =  lim 

t~+oo  t  t-^oo  t 

r  ^  1  ->  1- 

hm  ~  In  e  =  lim - 

t— ►OO  t  t-^OO  \j 


(2.38) 

(2.39) 


Obviously,  evaluating  equations  2.38  and  2.39  as  t  approaches  infinity  does  not  yield 
the  correct  exponents.  However,  the  derivation  of  the  extended  Lyapunov  exponents 
stipulated  that  $(to,to)  =  I.  Evaluating  and  $2  at  t  =  0  yields 


^i(to,to)  — 


^2(foi  to) 


0  -1 

1  0 


(2.40) 


(2.41) 


Thus,  defining  4&(to,  to)  =  I  is  absolutely  necessary  to  ensure  the  correct  propagation 
of  the  initial  conditions  at  t  =  0  to  the  coordinate-based  final  conditions  at  time  t. 
Therefore,  the  requirement  $(to,fco)  =  I  selects  a  unique  solution  to  $  =  A$,  and 
maintains  coordinate  continuity  from  t  =  0  to  time  t.  Additionally,  the  reader  should 
note  that  for  classical  Lyapunov  exponent  calculation,  any  solution  to  $  =  A$ 
suffices  to  determine  the  length  change  of  individual  vector  solutions. 

To  further  illustrate  the  necessity  of  coordinates,  consider  the  linear  constant 
coefficient  system  represented  by 


0  ijj 

—u)  0 


(2.42) 


Upon  inspection,  two  solutions  to  $  =  A$  are 


^i(t)  = 


coswt  sinwt 


—  sinwt  coswt 
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(2.43) 


$2(t)  = 


cos(a;t  +  ?/>)  sin(a;t  +  4') 
—  sin(u;t.  +  0)  cos(a;t  -f-  0) 


Both  matrices  have  column  solution  vectors  whose  norms  are  one  for  all  time.  Thus, 
the  classic  Lyapunov  exponents  are  zero  for  this  system,  and  the  eigenv'alues  of  A 
are  ±ia;  (purely  complex).  Therefore,  it  is  expected  that  the  extended  Lyapunov 
exponents  will  be  ±\u. 

Now,  to  calculate  the  extended  Lyapunov  exponents  for  the  system  defined  by 
$1,  the  first  thing  to  do  is  to  calculate  the  eigenvalues  of 


T  A( — 2  cos  tut)  T  1  —  0 


(2.44) 


2  cos  wt  ±  \/4cos2t4;t  —  4 

2 

=  cos<ut  ±  \/ cos^  tut  —  1 
=  cos  tut  ±  \/—  sin^tut 
=  cos  tut  ±  i  sin  tut 


(2.45) 


Using  these  values  the  extended  Lyapunov  Exponents  for  this  system  are  found 


to  be 


=  lim  -  ln(costut  ±  i  sintut) 

t->oo  t 

=  lim  -  Ine^'^'^’^^ 

t— oo  t 

lim  -(iftut)  =  ±itu 

t—oo  t 


(2.46) 


Hence.  produced  the  expected  values;  the  eigenvalues  of  the  A  matrix. 
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Similar  treatment  of  4^2  yields  a  dependence  on  i/’i 

iOi  =  lim  -(±/(a;t  +  I'y))  (2.47) 

t-<-oo  t 

If  0  =  0,  there  is  total  agreement  with  the  expected  values;  however,  for  0  =  — ut, 
the  extended  Lyapunov  exponents  equal  zero. 

Thus,  again  it  is  necessary  to  choose  a  coordinate  frame.  Evaluating  and 
4^2  at  t  =  0  produces 


^i(to,to)=^^  (2.48) 

(cos  0  sin  0  \ 

1  (2.49) 

—  sin  0  cos  0  y 


This  implies  that  in  order  to  ensure  correct  propagation  from  t  =  0  to  time  t,  the 
requirement  that  $(to,to)  =  I,  which  implies  0  =  0,  must  be  enforced. 

The  previous  examples  show  that  even  equation  1.9,  the  constant  coefficient 
system  solution,  is  a  coordinate- based  solution.  Additionally,  it  has  been  shown  that 
the  proposed  method  for  calculating  the  extended  Lyapunov  exponents  will  produce 
the  correct  system  roots,  both  the  real  and  imaginary  parts. 

These  examples  also  show  that  without  a  coordinate-based  description  of  mo¬ 
tion,  it  is  meaningless  to  ever  compute  an  oscillation  frequenc;, ,  the  imaginary  part 
of  the  appropriate  exponent,  for  even  a  constant  coefficient  system.  So,  there  should 
be  no  surprise  that  a  coordinate  frame  must  be  chosen  for  the  general  case. 


‘2.S  The  Limit  of  the  Mean 

There  are  times  when  the  extended  Lyapunov  exponent  functions  do  not 
necessarily  possess  limits  in  the  classical  sense,  as  time  goes  to  infinity.  This  section 
will  introduce  an  apparently  new  type  of  limit  called  the  limit  of  the  mean  (l.o.m.) 
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(11).  Tlu'ii,  it  will  He  argued  that,  in  the  sense  defined  by  the  l.o.in.,  the  limit  for 
those  extended  Lyapunov  exponent  functions  does  indeed  exist. 

Starting  out  with  the  definition  of  the  mean  value  of  a  function  on  the  interval 

(to.t): 

(f)(t)s7— !-t /'f(r)dT  (2.30) 

ft  —  to)  Jto 

Holding  the  lower  limit  to  constant,  the  limit  of  the  mean  will  be  defined  as 

l.o.m.f(t)  =  lim  (f)  (t)  =  liri  ^  -  /  f(r)dr  (2.51) 

t-oo  (t  -  to)  ito 

Apostol’s  limit  in  the  mean  (l.i.m.)  (1),  which  applies  to  the  convergence  of  a 
sequence  of  functions  to  a  limiting  function,  and  Cesaro’s  summability  for  an  infinite 
series,  are  conceptually  similar  to  the  limit  of  the  mean.  In  his  paper,  VViesel  (11) 
asserted  the  following  fundamental  result  and  proof: 
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Theorem:  For  any  inlegrable  f(t),  if  the  limit 


lim  f(t)  =  foo  (2.52) 

t— ►OO 

e.xists,  then  the  limit  of  the  mean  also  exists  and  equals  foo- 

If  linit— cof(t)  exists,  then  for  any  t  there  exists  a  T  such  that  when  t  >  T, 
|f(t)  ~  fco|  <  €■  Splitting  the  integral  of  2.47  yields 


l.o.m.f(t)  =  lim  /  f(r)dr+  lim  /  f(T)d7 

t— *oo  t— *oo  Ji  t— >00  J'Y 


(2.53) 


Assuming  f(t)  is  integrable,  the  first  integral  is  finite  and  constant  for  a  given  value 
T.  Thus,  its  limit  is  trivially  zero.  The  second  integral  is  bounded  by 


-  £)dr<  ^Xf(.)dr<  +  Od 


(t  ~  to) 


(t  —  to) 


(2.54) 


Performing  the  outer  integration  generates  the  following: 


(t-T) 
(t  —  to) 


Evaluating  2.55  as  t  — >  oo  yields 

(foo  -  e)  <  lim  77-^— T  /  f(r)dr  <  (foo  +  e) 
t-^oo  (t  —  to)  n 

Using  the  definition  of  the  limit  of  the  mean,  this  implies  that 

(foo  -  e)  <  l.o.m.f(t)  <  (foo  +  e) 

t— ►OO 

However,  since  this  is  true  for  any  e,  then 


l.o.m.f(t)  =  fo 

t— *•00 


(2.55) 

(2.56) 

(2.57) 

(2..58) 


Q.E.D. 
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Thus,  similar  to  the  concept  of  the  limit  in  the  mean,  the  limit  of  the  mean 
may  exist  when  no  classical  limit  exists.  To  illustrate  this,  it  is  easily  seen  that 

l.o.m.sint  =  0  (2.59) 

even  though  no  traditional  limit  exists.  The  limit  of  the  mean  extends  existence 
to  a  broader  class  of  functions,  but  not  an  unreasonably  broad  class.  For  example, 
the  limit  of  the  mean  does  not  exist  for  tsint,  nor  does  it  exist  for  any  monotonic 
unbounded  function. 

2.4  Extended  Lyapunov  Exponents  (Final  Form) 

In  section  2.1,  the  extended  Lyapunov  exponent  was  defined  by  equation  2.6. 
If  the  eigenvalues,  Aj(t),  of  the  fundamental  matrix  grow  exponentially  in  the  long 
term,  then  In  Aj  grows  no  faster  than  linearly  with  time.  Therefore,  a;,(t)  may  perhaps 
be  oscillatory,  but  it  will  have  no  long  term  tendency  to  grow.  Thus,  u;j(t),  in  all 
likelihood,  belongs  to  the  class  of  functions  for  which  the  limit  of  the  mean  exists. 

Hence.  Wiesel’s  (11)  final  conjecture:  The  classical  Lyapunov  exponent  is  the 
real  part  of 


w.(oo) 


l.o.m.u;,(<) 

i-*O0 


l.o.m. 

t— *00 


In  Ai(t) 

(t  -  to) 


In  A,(r) 


dr 


r 


(2.60) 
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where  the  limit  of  the  mean  must  be  calculated  over  extremal  pathways  through 


a  succession  of  eigenvalue  bifurcations,  and  the  imaginary  part  denotes  the  mean 
oscillation  frequency  with  respect  to  the  chosen  coordinate  frame.  Obviously,  the 
real  part  is  the  average  exponential  rate  of  increase  of  the  solution.  The  terms 
extremal  pathways  and  bifurcations  will  be  addressed  further  in  future  chapters. 

Validating  this  conjecture  is  the  main  thrust  of  this  thesis.  If  true,  it  could 
lead  to  a  unification  of  the  theory  of  linear  systems. 
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III.  Numerical  Algorithm  Development 


This  chapter  describes  the  algorithm  used  to  compute  the  extended  Lya¬ 
punov  exponents  for  second  order  systems,  and  the  associated  FORTRAN-coded 
programs/subroutines.  Generally  speaking,  a  predictor-corrector  numerical  integra¬ 
tion  scheme  (Haming)  is  used  to  perform  the  following  integrations; 


X 

=  F(X,t) 

(3.1) 

= 

(3.2) 

Wi 

=  -  j - ^^-^dr  (the  extended  Lyapunov  exponent) 

(3.3) 

The  main  program  reads  the  following  inputs:  initial  conditions,  system  parameters, 
the  maximum  time  of  integration,  and  the  number  of  steps,  and  it  sets  up  the  output 
files.  Next,  it  initializes  the  fundamental  matrix  $  to  the  identity  matrix,  and  the 
integrals  of  the  eigenvalues  to  zero.  The  “length”  of  the  vector  to  be  integrated 
is  then  defined  in  the  main  program,  and  the  timestep  is  calculated.  Then,  the 
numerical  integrator,  Haming,  is  initialized.  A  fourth  order  predictor-corrector  needs 
four  values  of  the  state  vector  to  perform  the  integration,  and  at  t  =  0,  there  is  only 
one  (the  initial  conditions).  Haming  uses  a  Picard  iteration  to  get  the  other  three 
values.  Provided  this  initialization  is  successful,  the  main  program  begins  the  full 
scale  integration  by  calling  Haming  again.  After  Haming  returns  each  value,  the  main 
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program  calls  the  subroutine  that  checks  for  eigenvalue  bifurcations  (SPLIT)  and 


the  subroutine  that  adjusts  the  phase  of  the  eigenvalues  (PHASE).  Finally,  the  main 
program  computes  the  actual  extended  Lyapunov  exponents.  All  of  the  required 
programs/subroutines  used  for  calculating  the  extended  Lyapunov  exponents  may 
be  found  in  Appendix  A,  and  are  summarized  in  Table  3.1. 


SUBROUTINE 

PURPOSE 

CALLED  BY 

Haming 

Numerical  integration 

Main  program 

SPLIT 

RHS 

Calculates  right  sides  of  equations  3.1  -  3.3 

Haming 

AMAT 

Calculates  the  elements  of  the  A  matrix 

RHS 

CALEIG 

Calculates  complex  logs  of  eigenvalues  of  $ 

RHS 

SPLIT 

Checks  for  bifurcations  and  if  one 
has  occurred,  restarts  Haming 

Main  program 

SPLIT2 

Determines  bifurcation  time  and  type 

SPLIT 

Extrapolates  log  eigenvalues 
into  and  out  of  bifurcation 

Swaps  root  labels  if  necessary 

PHASE 

Fixes  phase  on  current  eigenvalues  if 
necessary 

Main  program 

SPLIT2 

CHECKS 

Performs  reasonableness  checks  on 
eigenvalues  and  averaged  integrals 

SPLIT 

Table  3.1.  Subroutines  for  Calculating  the  Extended  Lyapunov  Exponents 


The  first  step  in  he  algorithm  is  to  set  up  the  equations  to  be  integrated  3.1 
to  3.3.  The  next  three  sections  discuss  this  setup  in  detail  for  second  order  systems. 
Section  4  discusses  the  development  of  the  subroutine  to  calculate  the  eigenvalues  of 
the  ^  matrix,  and  the  last  .section  addresses  potential  problems. 
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3.1  Deriving  the  General  Form 


Starting  with  a  second-order  ordinary  differential  ecjuation,  one  must  transform 
it  into  the  form  X  =  F(X,t).  For  example,  let  the  following  represent  the  original 
differential  equation  where  a,  b  may  be  functions  of  x  and/or  t. 


X  -(-  ax  -f-  bx  =  0 


(3.4) 


Now,  let 


Xi  =  X 

=»  Xi  =  X 

and 

X2  =  Xi 

^X2  =  X  =  (— ax2  -  bxi) 

This  yields, 

Xl 

X2 
or, 

Xl  =  X2 

X2  =  — bxi  —  ax2 


' 

0  1 

> 

Xl 

— b  —a 

X2 

=  F(X,t) 


(3.5) 
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3.2  Calculating  the  Right-Hand  Sides  (RHS) 


The  right-hand  side  of  3.6  is  calculated  in  the  subroutine  RHS  (see  Appendix 
A),  which  is  called  by  flaming.  Additionally,  RHS  calculates  the  terms  in  integral 
3.2.  The  calculation  of  the  A  matrix  is  discussed  in  the  following  section,  and 
is  performed  in  the  subroutine  AM  AT;  however,  RHS  calls  AM  AT,  and  then  the 
product  of  the  matrices  A  and  $  is  performed  in  RHS  and  thus,  is  outlined  here. 


Thus, 


$  =  A^ 


.. 

d>ii 

4>vi 

_ 

an  ai2 

4>12 

<f>2l 

4>22 

a2i  a22 

<t>2\ 

(f)22 

(3.6) 

(3.7) 


<?ii  =  (aii(/>ii  +  a.i2<?!>2] )  (3-8) 

—  (^11(^12  +  ai2<^22) 
d>21  =  (a21</’ll  +  ft22^2l) 

O22  —  (321^12  +  <^22^22) 


Kfiuations  3.8  are  evaluated  in  the  .subroutine  RHS. 

rhe  last  thing  that  RHS  does  is  to  take  the  instantaneous  time-average  of  the 
logarithms  of  the  eigenvalues  by  dividing  In  A,  by  the  instantaneoiis  t.  The  actual 
calculation  of  the  eigenvalues  A,  of  the  fundamental  matrix  $  is  jK'rformed  in  the 


subroutine  CALEIG  which  is  described  in  section  3.4,  and  which  is  called  by  RHS. 


3.3  Finding  the  .Associated  Linear  System  (AM.AT) 

In  order  to  perform  the  integration  of  3.2,  the  A  matrix  must  first  be  computed. 
This  relatively  simple  calculation  is  performed  in  the  subroutine  AMAT  according 
to  the  following  development. 

Recall  A(t)  was  defined  as  a  Jacobian  matrix;  the  matrix  of  partial  derivatives 
of  the  system  differential  equations  with  respect  to  the  variables  of  the  problem. 
Therefore, 


d%7 

3x7 

0 

1 

A(t)  = 

d%\ 

3X2 

= 

9(— bxt  —3x2) 
axi 

9(-bxi-ax2) 

3X2 

>< 

1 

_ _ 1 

(-a . 

(3.9) 


The  subroutine  RHS  calls  AMAT  just  prior  to  calculating  the  product  A$. 


3.4  Eigenvalues  of  the  Fundamental  Matrix  (CALEIG) 

The  eigenvalues  of  a  2  x  2  fundamental  matrix  are  easily  calculated  from  its 
characteristic  equation,  which  is  defined  by: 


det  -  AI} 


n2 


d>2\  <t>22  ~  ^ 


=  0 


(3.10) 
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Evaluating  the  determinant  yields 


[((/'ll  ~  —  A)  —  </>i2</'2l]  =  0 


(3.11) 


A^  -f  A(  — ^11  +  </'22)  +  (<^ll</'2'2  ~  </'l2^2l)  —  0 


(3.12) 


Solving  the  quadratic  produces  the  following  equation  for  A: 


((/’ll  +  ^22)  i  J{  —  4>\\  ~  4>22Y  ~  4(^11</»22  —  4>\2<t>2\) 


(3.13) 


In  addition  to  calculating  3.13,  the  subroutine  CALEIG  computes  the  complex 
logarithm  of  the  eigenvalues  and  returns  them  to  RHS  separated  into  their  real  and 
imaginary  parts.  To  implement  this  separation,  one  needs  to  evaluate  the  natural 
logarithm  of  a  complex  number.  Let  the  following  equation  define  an  arbitrary 
complex  number: 


z  =  Q  -f 


—  re'^  (polar  form) 


where. 


r  =  l^l  =  +  iP 


(3.14) 
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and 


9  =  arg  z  =  arctan 

Taking  the  logarithm  of  z  in  its  polar  form  yields 

ln(re‘^)  =  lnr  +  lne‘^  (3.15) 

=  In  r  +  iO 

Thus,  the  real  part  of  In  2  is  In  r,  and  the  imaginary  part  is  0. 

These  are  the  values  returned  to  RHS  by  CALEIG.  RHS  then  divides  the 
natural  logarithms  of  the  eigenvalues  by  t  to  yield  the  instantaneous  rates  of  change. 

3.5  Potential  Problems 

As  was  alluded  to  in  the  introductory  part  of  this  chapter,  there  are  situations 
when  the  calculation  of  the  eigenvalues  of  the  fundamental  matrix  runs  into  problems. 
One  of  these  problems  is  the  fact  that  the  eigenvalues  are  growing  exponentially,  and 
thus  can  quickly  span  several  orders  of  magnitude.  As  this  occurs,  floating  point 
overflow  or  underflow  can  result  in  the  calculations. 

The  second  problem  of  significance  that  can,  and  does,  arise  is  eigenvalue  bifur¬ 
cation:  two  real  eigenvalues  merging  into  two  complex,  or  two  complex  merging  into 
two  reals.  Figure  3.1  shows  a  typical  bifurcation.  Typically,  there  will  be  an  end¬ 
less  succession  of  these  eigenvalue  bifurcations,  with  the  time  between  bifurcations 
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decreasing  rapidly. 


Re[lii  X,]  - 


25  3  35  4 

Time 


Figure  .3.1.  Sample  Bifurcation 

First,  the  software  must  be  able  to  detect  that  an  eigenvalue  bifurcation  has 
occurred.  The  bifurcation  point  is  relatively  straightforward  to  find  for  second-order 
.systems,  because  the  discriminant  of  the  characteristic  equation  will  change  sign 
when  a  bifurcation  occurs.  Thus,  comparison  of  the  last  discriminant  calculated 
with  the  present  one  will  detect  whether  or  not  a  bifurcation  has  occurred.  The 
subroutine  SPLIT  performs  this  comparison. 

If  a  bifurcation  has  occurred,  then  the  softw'are  must  integrate  through  the 
bifurcation,  and  label/identify  the  eigenvalues  as  they  merge.  This  is  performed 
in  the  subroutine  SPLIT2,  which  is  called  from  SPLIT.  Figure  3.1  illustrates  that 
the  slopes  of  the  eigenvalues  are  infinite  at  the  time  of  bifurcation.  Defining  a 
new  time  variable  .s  =  \t  —  ,  where  tb  is  the  bifurcation  time,  eliminates  this 

difficulty.  The  bifurcation  time  is  calculated  in  SPLIT2  by  noting  that  time  is  a 
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locally  ciuadratic  function  of  the  difference  between  the  two  merging  eigenvalues  (real 
or  imaginary  parts  as  appropriate).  Thus,  the  time  of  bifurcation  is  determined  by 
inverse  interpolation  back  to  the  time  of  the  bifurcation.  With  the  bifurcation  time 
known,  the  merging  eigenvalues  are  determined  by  extrapolating  in  the  s  variable  to 
s  =  0.  (Polynomial  extrapolation  in  time  is  impossible  because  of  the  infinite  slopes 
seen  in  Figure  .3.1.)  Further  extrapolation,  .slightly  away  from  .s  =  0  (still  using  the 
s  variable),  enables  the  numerical  integration  to  be  restarted.  This  entire  process 
introduces  only  minimal  errors,  and  occurs  within  one  normal  integration  timestep. 

It  is  essential  that  the  software  keep  track  of  which  eigenvalue  is  which  during  a 
bifurcation  so  that  characterization  of  their  long  term  behavior  is  not  compromised. 
Thus,  SPLIT2  determines  the  bifurcation  type  (real-to-imaginary,  or  imaginary-to- 
real),  and  tracks  the  eigenvalues  according  to  this  type.  If  the  type  is  determined  to 
be  a  pair  of  reals  going  to  complex,  SPLIT2  insures  that  the  root  with  the  positive 
imaginary  part  emerges  from  the  bifurcation  with  a  positive  slope.  Since  the  two 
merging  eigenvalues  lose  their  identity  in  the  bifurcation,  SPLIT2  is  free  to  swap  root 
labels  if  necessary.  If  the  type  is  a  pair  of  complex  values  merging  to  real,  SPLIT2 
calls  the  subroutine  PHASE  to  reset  the  27r  multipliers  of  the  eigenvalues.  It  is  only 
after  this  labelling,  or  tracking,  has  occurred  that  SPLIT2  extrapolates  the  average 
eigenvalues  over  the  bifurcation. 

The  PHASE  subroutine  mentioned  previously  is  also  called  from  the  main 
program  to  update  the  phase  of  the  eigenvalues  as  they  are  propagated  forward. 


I'his  ensures  that  the  imaginary  part  of  In  (A,)  will  have  a  net  increase. 

An  additional  check  on  the  eigenvalues  for  second  order  systems  is  the  conju¬ 
gality  of  the  eigenvalues.  In  the  subroutine  CHMCKS,  the  imaginary  parts  of  the 
eigenvalues  and  their  averaged  integrals  are  checked  for  pairedness.  Specifically,  the 
absolute  value  of  the  sum  of  the  imaginary  parts  should  be  zero,  within  numerical 
accuracy. 

The  algorithms  implemented  to  address  these  problems  will  be  successful,  if 
the  growth  of  the  eigenvalues  stays  reasonable  during  the  time  of  integration,  and 
the  bifurcations  do  not  occur  back-to-back  during  a  single  timeste]).  If  the  growth  is 
too  rapid,  accuracy  is  lost,  and  if  back-to-back  l)ifurcations  occur,  (he  bifurcations 
will  not  even  be  detected.  In  either  of  these  cases,  characterization  of  the  long  term 
behavior  will  be  compromised. 


;mo 


IV.  Application  and  Validation 


Application  of  the  technique,  developed  in  previous  chapters,  to  the  van  der 
Pol  equation  is  offered  here.  The  first  sections  contain  the  .setup  and  anal,ysis  of 
different  trial  cases.  A  Fourier  analysis,  used  for  validation  of  the  imaginary  parts 
of  the  extended  Lyapunov  exponents,  is  pre.sented  in  the  last  section. 

/,.I  \  'an  der  Pol's  equation 

Van  der  Pol's  equation  is  a  second-order  differential  equation  with  a  linear 
restoring  force  and  nonlinear  damping,  that  originally  arose  as  an  idealization  of  a 
self-excited,  or  spontaneously  oscillating,  valve  circuit(5).  Per  .Jordan  and  Smith  (.5), 
the  van  der  Pol  equation  may  be  written  a.s 

X  +  e(x''^  -  l)x  +  X  =  0  (4.1) 

Figures  4.1  and  4.2  contain  the  phase  portraits  for  the  van  der  Pol  equation's 
two  stable  configurations.  For  c  <  0.  the  van  der  Pol  equation  has  a  stable  equi¬ 
librium  point  at  the  origin.  As  seen  in  Figure  4.1,  all  nearby  trajectories  spiral  in 
toward  the  origin.  Figure  4.2  illustrates  that,  with  an  f  >  0,  the  van  der  Pol  equation 
exhibits  a  limit  cycle  behavior,  and  is  considered  the  cla.ssic  paradigm  for  self-excited 


oscillations(7). 


4.1.1  Sehip  First,  this  must  be  put  into  the  general  form  X  =  F(X,  t).  To 


do  this,  let 


Xi  =  X 

X2  =  Xi 


which  yields, 

Xi  =  X2 

and  (4.2) 

X2  =  -e(x^  -  1)X2  -  Xi 


These  are  the  equations  required  for  the  subroutine  RHS. 

Next,  the  equations  for  calculating  the  A  matrix  must  be  setup.  Recall, 


A  = 


dF 


(4.3) 


Evaluating  equation  4.3  produces  the  following  expressions  for  the  elements  of  A: 


an  = 

ai2  = 

^21  = 


dx2 


dx 


dx2 

dx2 


d[-({x\  -  1)X2  -  Xi 
dX] 


=  — 2cXiX2  —  1 


(4.4) 
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^i22 


g[-e(x^  -  1)X;;  -  Xi] 

dx2 


=  -e(x'i  -  1) 


These  are  the  equations  that  will  be  coded  into  the  subroutine  AMAT. 


4.1.2  Case  1  (The  Constant  Coefficient  Case)  Linearizing  about  the  equi¬ 
librium  point,  which  is  located  at  the  origin  for  the  van  der  Pol  equation,  yields  an 
A  matrix  of  constant  coefficients  for  the  equation  $  =  A$,  and  the  eigenvalues  of 
A  should  correspond  to  the  calculated  extended  Lyapunov  exponents.  Then  the  A 
matrix  at  t  =  0  is  equal  to  the  A  matrix  at  all  time: 


A(t)  =  A(to)  = 


0  1 
-1  e 


Solving  for  the  eigenvalues  of  A 


det  {A  —  AI} 


det 


0-  A 


1 


—1  c  —  A 
A^  -  eA  -f  1  =  0 


=  0 


(4.5) 


(4.6) 

(4.7) 


Solving  equation  4.6  yields 

A  - 

2 


(4.8) 
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From  equation  4.8,  it  is  abvious  that  choosing  |cj  small  (i.e.  |e|  <  2)  will  ensure  an 
imaginary  part  of  the  eigenvalue.  For  this  analysis,  e  =  —  1  will  be  used.  Therefore, 


A 


-1  ±  \/n^ 


-0.5 

2 


-0.5  ±  (0.866)i 


(4.9) 


Figure  4.3  contains  the  phase  portrait.  Note  that  the  equilibrium  point  is  stable,  and 
the  trajectory  is  a  spiral  towards  the  point  (0,0).  Jordan  and  Smith’s  (5)  treatment 
of  the  van  der  Pol  equation  contains  a  proof  that  the  requirement  for  stability  of  the 
equilibrium  point  is  a  e  <  0.  Since  e  =  —  1  was  chosen,  the  phase  portrait  shows  that 
the  system  behaves  in  the  expected  fashion. 


Figure  4.3.  Van  der  Pol  Equation  Phase  Portrait  (Stable  Equilibrium  Point  Case) 
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Figures  4.4  and  4.5  show  the  calculated  real  parts  of  the  In  A,  versus  time.  As 
expected,  the  logarithms  of  the  eigenvalues  are  changing  linearly  with  time,  with  a 


slope  of  -0.5. 


Re[ln 


Time 


Figure  4.4.  Real  Part  of  In  Ai  versus  Time 
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Figure  4.5.  Real  Part  of  In  A2  versus  Time 
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F'igure  4.6  shows  the  calculated  imaginary  parts  of  In  A,  versus  time.  The 
imaginary  parts  are  changing  linearly  with  time  (sloped  0.866).  until  t  .3.6.  At 
this  time,  the  In  A,  undergo  back-to-back  bifurcations  that  cannot  be  detected  by  the 
software. 


1  m[  1  n  X|  J  4 


i 

1 


0  12  3  4 


Tl  me 


Figure  4.6.  Imaginary  Parts  of  In  A,  versus  Time 
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F'igure  4.7  is  a  plot  of  the  characteristic  ecjuation’s  discriminant  versus  time. 


This  plot  illustrates  why  there  is  a  problem  at  t  ss  3.6.  For,  at  this  time,  the  back-to- 
back  bifurcations  occur  so  close  together  that  the  discriminant  never  appears  to  even 
change  sign.  Thus,  the  software  never  detects  that  a  bifucation  has  occurred.  This 
situation  did  not  improve  even  when  the  program  was  run  at  the  smallest  timestep 
that  was  possible. 
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Figure  4.7.  Discriminant  of  the  Characteristic  Equation 
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riius.  luimorical  iiiaccurarios  force  tlie  integration  to  be  limited  to  before  the 
principal  bifurcation.  It  is  obvious  though  from  figures  1.7  and  1.8.  a  i)lot  of  one  of 
the  variabh's  versus  time,  that  the  solution  quickly  converges  to  a  s7c«dy-.s7«/r  vaiur 
in  this  ca.se.  Also,  as  the  following  figures  will  indicate,  e.xcellent  results  were  still 
obtainable. 


Time 


Figure  4.8.  X)  versus  I’ime 


Figures  4.9  and  4.10  contain  plots  of  the  real  parts  of  the  calculated  extended 


Lyapunov  exponents.  Clearly,  even  with  the  bifurcation  problem  that  was  encoun¬ 
tered,  they  are  conv^erging  to  the  values  obtained  for  the  real  parts  of  the  eigenvalues 
of  A:  -0.5  ( Reference  equation  4.9). 
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Figure  4.9.  Real  Part  of  Extended  Lyapunov  Exponent  1 
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Figure  4.10.  Real  Part  of  Extended  Lyapunov  Exponent  2 
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Figure  1.11.  the  plot  of  the  imaginary  parts  of  the  extended  Lyapunov  ex¬ 
ponents,  again  shows  that  the  technique  has  yielded  the  correct  values.  Clearly,  as 
t  — >  c»,  the  imaginary  i)arts  are  converging  to  some  number  close  to  ±0.86.  This  is 
tlie  same  value  that  was  calcidatenl  for  the*  imaginary  parts  of  the  eigenvalues  of  the 
constant  coefficient  A  matrix. 
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Figure  4.1 1.  Imaginary  Part.s  of  Exlfiidfd  Lyapunov  Exponents 


It  should  be  noted  here  that  because  the  data  was  suspect  after  the  principal 


bifurcation,  it  was  not  included  for  the  extended  Lyapunov  exponent  plots.  However. 


the  plots  flo  conclusively  show  that  the  technique  has  yiekh'd  the  correct  values  for 


the  system  characteristic  exponents. 


Case  2  (The  Periodic  Coefficient  Case)  Linearizing  about  the  limit 
cycle  yields  an  A  matrix  of  periodic  coefficients  for  the  ecjuation  $  =  A$.  Thus, 
the  Poincare  exponents  for  this  case  should  correspond  to  the  extended  Lyapunov 
exponents  calculated  using  the  technique  outlined  in  the  thesis.  Per  Jordan  and 
Smith  (5),  the  limit  cycle,  for  c  =  1,  can  be  defined  by  the  following  initial  conditions: 


Xi(to)  ^  2 

X2(to)  ~  0 

The  angular  frequency  of  the  limit  cycle(5)  is 

27r 

LJ  =  - 

T 


(4.10) 

(4.11) 


(4.12) 


where  r  is  the  period  of  the  limit  cycle. 

According  to  Floquet  theory,  the  Poincare  exponents  are  defined  as 

p,  =  -  In  A,(r)  ±  (4.13) 

r  r 


where  r  is  the  period  of  the  trajectory  and  A,(r)  are  the  eigenvalues  of  the  fun¬ 
damental  matrix  #  at  t  =  r.  Thus,  in  order  to  calculate  the  Poincare  exponents, 
the  period  of  the  limit  cycle  must  be  determined,  and  then  the  eigenvalue  of  the 
fundamental  matrix  at  the  end  of  one  period  in  time  must  be  calculated. 
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From  Figure  4.12,  which  is  the  plot  of  Xi  versus  time,  it  can  be  seen  that  the 


period  is  equal  to  approximately  6.67. 
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Figure  4.12.  Xi  versus  Time 
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Integration  of  the  fundamental  matrix  $  through  this  time  gave  the  following 


values  for  $(r): 


$(r) 


0.01998446583557  2.685067523504 

0.006971779004636  0.9788204972462 


(4.14) 


A  restatement  of  equation  3.13,  which  is  used  to  calculate  the  eigenvalues  of 
the  fundamental  matrix  follows: 

{4>\\  4>22)  ^  —  —  <t>22Y  —  ~ 

A  =  - ^ -  (4.15) 


Evaluating  the  above  using  the  elements  of  the  $  matrix  from  equation  4.14  yields 
the  following  values  for  Aj(r) 


Ai(r)  =  0.997962 

(4.16) 

A2(r)  =  .000843226 

(4.17) 

(4.18) 

The  Poincare  exponents  are  then  calculated  using  these  values  for  Ai(r)  in  equation 
4.13,  and  they  are  found  to  be 


Pi 


P2 


In  A2 


-1.05 


(4.19) 

(4.20) 
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rh('  leader  should  note  that  tlie  db?27r/r  term  of  equation  4. Id  has  been  dropped. 
This  is  standard  practice  in  Floquet  theory  because  this  term  is  the  uncertainty  that 
arises  due  to  the  limits  on  the  natural  logarithm  function. 
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Figure  4.13  contains  the  phase  portrait  which  was  plotted  with  the  system 


elements  calculated  using  the  technique  outlined  in  the  previous  chapter,  for  the 
following  initial  conditions  and  system  parameters: 


xi  =  2.0 

X2  =  0.0 

e  =  +1 


Obviously,  the  trajectory  is  the  limit  cycle  given  in  Jordan  and  Smith  (5). 


Figure  4.14  contains  the  plot  of  the  real  parts  of  In  A,  versus  time.  As  expected, 
the  logarithms  are  changing  linearly  with  time,  and  engage  in  a  continual  series  of 
bifurcations.  These  bifurcations  begin  to  cause  problems  at  about  t  =  6  because  of 
the  small  amount  of  time  between  bifurcations.  Additionally,  it  can  be  seen  that  the 
eigenvalues  themselves  are  now  spread  over  at  least  five  orders  of  magnitude,  with 
the  smallest  being  very  close  to  zero.  Thus,  accuracy  is  compromised  at  the  end  of 
the  integration  time. 
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Figure  4.14.  Real  Parts  of  In  Aj  versus  Time 
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Figure  4.15  shows  the  calculated  imaginary  parts  of  the  InA^  versus  time.  The 
imaginary  parts  show  a  characteristic  (11)  “stairstep”  behavior  with  the  phase  angle 
continually  increasing,  on  the  average.  Recall,  the  identity  of  each  root  is  chosen  at 
the  time  of  bifurcation  to  ensure  that  this  occurs.  This  is  absolutely  necessary  if  the 
long  term  behavior  of  the  system  is  to  be  characterized. 


Figure  4.15.  Imaginary  Parts  of  In  A,  versus  Time 
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For  comparison’s  sake,  Figures  4.16  and  4.17  are  plots  of  the  discriminant  of 
the  characteristic  equation  for  this  case.  It  is  obvious  from  4.16  that  numerous 
bifurcations  are  ocurring  and  that  after  the  principal  bifucation,  they  are  ocurring 
closer  in  time. 
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Time 


Figure  4.16.  Discriminant  of  the  Characteristic  Equation 
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Figure  4.17  shows  that,  especially  for  the  principal  biuircation,  the  discriminant 
changes  sign,  and  remains  negative  for  a  time.  This  implies  that  there  is  not  a  back- 
to-back  bifurcation  at  the  principal  bifurcation  time.  Thus,  the  software  detected  the 
bifurcation  and  successfully  tracked  the  eigenvalues  through  the  bifurcation,  unlike 
the  bifurcations  in  the  equilibrium  point  case. 


Figure  4.17.  Discriminant  of  the  Characteristic  Equation 
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The  real  parts  of  the  extended  Lyapunov  exponents  versus  inverse  time  are 
plotted  in  Figures  2.3  and  2.4.  Extrapolating  out  to  where  inverse  time  is  zero  (i.e.. 
t  ex),  the  real  parts  are  obviously  converging  to  the  real  parts  of  the  Poincare 
exponents:  0  and  -1.05. 
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Figure  4.18.  Real  Part  of  Extended  Lyapunov  Exponent  1 
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l'igur('  4.10.  Real  Part  of  Extendid  Lyapunov  ILxponent  2 


The  imaginary  parts  of  the  extended  Lyapunov  exponents  are  plotted  in  Fig¬ 
ure  4.20  below.  Although  calculation  of  the  Poincare  exponents  did  not  yield  an 
imaginary  part,  Figure  4.20  confirms  that  there  is  an  imaginary  part.  Though  the 
extended  Lyapunov  exponents  start  out  at  zero,  as  inverse  time  approaches  zero 
(t  ^  oo).  the  imaginary  parts  are  converging  to  u;,  ^  ±0.9,  which  is  approximately 
equal  to  the  '‘uncertainty’'  term  in  the  equation  for  calculating  the  Poincare  expo¬ 
nents:  ±ii^  =  ±^  =  0.94. 
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Figure  4.20.  Imaginary  Parts  of  Extended  Lyapunov  Exponents 

Thus,  although  standard  Floquet  theory  does  not  yield  an  imaginary  part,  the 
technique  outlined  in  this  paper  confirms  such  a  part  exists.  Further  validation  of 
existence,  and  correct  ness  of  the  valvies  obtained  for  the  imaginary  parts  is  present  ed 
in  the  next  .section. 
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4-2  Validation  of  the  Imaginary  Part  of  the  Extended  Lyapunov  Exponent 

To  validate  the  imaginary  part  of  the  extended  Lyapunov  exponent,  EFT  (fast 
Fourier  transform)  techniques(  12)  will  be  used  to  obtain  the  PSD  (power  spectral 
density)  plot  for  the  system  in  each  of  the  trial  cases.  Appendix  B  contains  the 
computer  code  used  for  this  analysis.  Most  often,  waveforms/data  are  measured 
in  the  time  domain.  Time-based  measurements  have  historical  precedence.  For 
example,  Galileo  reportedly  used  his  own  pulse  as  a  timepiece  in  making  his  original 
pendulum  observations.  Llnquestionably,  time-based  measurements  are  the  most 
familiar  data  format,  but  time  histories  tell  only  one  side  of  the  story.  The  type  of 
power  vs  frequency  spectrum  that  a  system  possesses  will  classify  its  response.  For 
instance,  if  the  spectrum  contains  n-number  of  spikes  at  discrete  frequencies,  it  is 
considered  a  periodic  response  of  period  n.  and  if  the  spectrum  is  continuous,  with 
no  discernible  spikes,  it  is  classified  chaotic.  Additionally,  if  there  are  spikes  in  the 
spectrum,  these  spikes  will  occur  at  the  driving  frequencies  of  the  system. 

4.2.1  .Application  of  EFT  Technique.^  For  each  of  the  trial  cases,  the  trajec¬ 
tories  were  numerically  integrated  using  4.2,  and  a  single  displaced  trajectory  on  the 
coordinatized  tangent  space  was  also  integrated  using 

<^X(t)  =  A(t)6X(t)  (4.21) 
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At  frequent,  equally  spaced  time  increments,  the  displacement  vector  6X(t)  was  nor¬ 
malized,  coordinate  values  were  saved,  and  the  integration  was  restarted.  Given  a 
random  initial  condition  for  ^X(to),  it  is  expected  that  the  trajectory  will  strongly 
converge  towards  the  extended  Lyapunov  exponent  with  the  largest  real  part.  Nor¬ 
malizing  the  vector  6X(t)  discards  the  exponential  growth  information,  and  retains 
only  information  on  how  the  vector  rotates  in  space  (i.e.,  it  retains  the  information 
on  the  imaginary  part).  Then,  individual  coordinate  histories  were  processed  with 
the  fast  Fourier  transform  to  extract  frequency  information.  The  subsequent  one 
sided  power  spectral  density  will  show  whether  the  relative  motion  has  a  discrete  or 
continuous  frequency  spectrum,  and  if  there  are  discrete  spikes,  these  will  occur  at 
the  fundamental  frequency  of  the  system. 
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4-2.2  Validation  of  the  imn  Her  Pol  Equilibrium  Point  Case  Application  of 
this  technique  to  the  van  der  Pol  equilibrium  point  case  produced  Figure  4.21. 
Clearly,  the  spectrum  is  discrete,  with  a  large  spike  at  angular  frequency  approxi¬ 
mately  0.9  and  a  much  smaller  peak  at  three  times  this  value. 
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h’igure4.21.  Power  Spectral  Density  of  a  Relative  Trajectory  versus  Angular 
Frequency 
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Figure  4.22  is  an  expanded  view  of  the  portion  of  the  spectrum  close  to  the 
largest  spike,  and  shows  that  the  large  spike  actually  occurs  at  something  less  than 
0.9,  about  0.86. 
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Figure  4.22.  Power  .Spectral  Density  of  a  Heiative  Trajectory  versus  Angular 
Frequency 


riius.  the  spike  occurs  at  exactly  the  value  ol)tained  for  the  imaginary  part  of 
t  he  extended  Lyapunov  exponent  (^5;  0.86),  and  the  imaginary  part  of  the  eigenvalues 
of  the  A  matrix  (0.86). 
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4.2.3  Validation  of  the  van  der  Pol  Limit  Cycle  Case  Applicat  ion  of  the  fast 
Fourier  transform  technique  to  the  van  der  Pol  limit  cycle  case  produced  Figure  4.23. 
Again  the  spectrum  is  discrete  with  a  large  spike  at  angular  frequency  approximately 
0.9  and  a  much  smaller  peak  at  three  times  this  value. 
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Figure  4.23.  Power  Spectral  Density  of  a  Relative  Trajectory  versus  Angular 
Frequency 
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Figure  4.24  is  an  expanded  view  of  the  portion  of  the  spectrum  close  to  the 
largest  spike,  and  shows  that  the  large  spike  actually  occurs  at  something  greater 
than  0.9,  about  0.94. 
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P'igure  4.24.  Power  Spectral  Density  of  a  Relative  Trajectory  versus  Angular 
Frequency 

.Again,  the  value  obtained  for  the  imaginary  part  of  the  extended  Lyapunov 
exponent  is  validated.  Additionally,  this  plot  serves  as  further  confirmation  that  an 
imaginary  part  of  the  Poincare  exponents  exists,  and  is  equal  to 


V.  Conclusions  and  Recommendations 


An  explicit  numerical  technique  has  been  developed  to  calculate  the  charac¬ 
teristic  exponents  {extended  Lyapunov  exponents)  for  any  second  order  system.  The 
tact  that  this  technique  can  be  applied  to  any  system  unifies  the  theory  of  second- 
order  linear  systems.  A  summary  of  conclusions,  and  recommendations  for  further 
investigation  are  presented  in  the  following  sections. 

5.1  Conclusions 

The  extended  Lyapunov  exponents  calculated  for  the  van  der  Pol  equation 
stable  equilibrium  point  case  were  in  excellent  agreement  with  the  theory:  they  ex¬ 
actly  equalled  the  eigenvalues  of  the  constant  coefficient  A  matrix.  Additionally,  the 
discrete  spike  in  the  PSD  at  the  value  of  the  imaginary  part  of  the  extended  Lya¬ 
punov  exponent  for  this  system  clearly  validated  the  imaginary  part  of  the  exponent 
as  being  the  system’s  fundamental  frequency. 

’[’he  value  of  the  extended  Lyapunov  exponent  for  the  periodic  case  (the  van 
der  Pol  equation  limit  cycle)  went  beyond  simple  agreement,  it  confirmed  the  exis¬ 
tence  of  an  imaginary  part  to  the  Poincare  exponent,  and  determined  its  value,  a 
question  that  Floquet  theory  left  unanswered  in  this  case.  Also,  the  clear  spike  in  the 
PSD  at  the  same  value  as  the  imaginary  part  of  the  extended  Lyapunov  exponent 


further  validated  its  value  and  existence. 


These  results  suggest  that  the  potential  for  unifying  linear  system  theory  in 
general  is  very  great.  The  concept  that  a  system’s  characteristic  exponents  contain 
both  a  real  and  imaginary  part,  whether  the  system  is  constant,  periodic,  multi¬ 
ply  periodic,  or  aperiodic,  is  not  only  intuitively  correct,  it  has  been  shown  to  be 
numerically  correct  for  the  trial  cases  presented  here. 

5.2  Recommendations 

Recommendations  for  further  study  of  the  extended  Lyapunov  exponent  and 
applications  of  this  technique  are: 

1)  Apply  the  technique  to  the  van  der  Pol  system  with  a  periodic  forcing 
function. 

2)  Apply  the  technique  to  other  classic  second-order  systems  (c.g..  the  pendu¬ 
lum). 

3)  Develop  the  code  and  apply  the  technique  to  some  classic  higher  order 
systems. 

The  major  difficulty  anticipated  with  these  recommendations  has  to  do  with  the 
calculation  of  the  eigenvalues  of  the  fundamental  matrix  The  timespan  over  which 
this  technique  is  applicable  is  currently  limited  by  pairs  of  eigenvalue  bifurcations 
that  are  separated  by  short,  exponentially  decreasing  time  intervals.  In  his  original 
paper  on  ext<'nded  Lya|)unov  exponents  (II).  Wiesel  developed  an  algorithm  for 
pro[)agating  the  eigenvalue  and  eigenvector  matrices  numerically.  In  essence,  it  was 
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shown  that  the  eigenvectors  and  eigenvalues  obey  differential  equations  of  their  own, 
thereby  eliminating  the  need  to  propagate  the  fundamental  matrix  itself.  Use  of  this 
algorithm  could  possibly  reduce  the  eigenvalue  bifurcation  difficulty.  Additionally, 
for  higher  order  systems,  direct  calculation  of  the  eigenvalues  of  the  fundamental 
matrix  $  becomes  much  messier  than  the  second  order  systems;  thus,  numerically 
integrating  the  eigenvalue  and  eigenvector  differential  equations  themselves  will  most 
likely  be  more  efficient  and  accurate  in  the  application  of  the  technique  to  higher 


order  systems. 


Appendix  A.  Computer  Programs  for  Calculating  the  Extended 


Lyapunov  Exponents 
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c 

c 

c 

c 

c 

c 


c 

c 

c 


program  vplmain 

2D  simple  van  der  Pol  system 
calculate  eigenvalues  from  phi 
and  time  averaged  eigenvalues 

common  /ham/  t,x(10,4) ,f (10,4) ,err(10) ,nn,hh,mode,e 
double  precision  t,x,f ,err,hh,e 

common  /eval/  eig(2,4) ,npi(2) ,discl(4) ,dtbif ,iswap(2) 
complex* 16  eig 

double  precision  disci, dtbif 

common  /debug/  idebug 

double  precision  tmax , omegr , omegi 

double  precision  a(2,2) ,phi(2,2) ,work(40) ,wpp(2) ,vmag 
double  precision  wqq(2) ,wrr(2) ,wss(2) 
complex* 16  wq,wr,ws 
complex*16  w(2) ,vec(2,2) ,wp 

equivalenceC  wp,wpp) , (wq,wqq) , (wr,wrr) , (ws,wss) 

input  parameters,  Ic’s 

read  (*,*)  t 
read  (*,*)  e 


c 

c  open  output  file 
c 

open ( 13 , FILE= ’ run . out ’ , STATUS= ’ UNKNOWN ' ) 

C  PRODUCTION  CODE 

OPEN (11, FILE= ’ lyap . re ’ , STATUS= ’UNKNOWN ’ ) 
OPEN ( 12 , FILE= ’ lyap . im ’ , STATUS= ’UNKNOWN ’ ) 

C  END  PRODUCTION  CODE 
C  CHECKOUT  CODE 

OPEN ( 1 , FILE= ’ Ineigl . re ’ , STATUS= ’ UNKNOWN ’ ) 
OPEN ( 2 , FILE= ’ Ine ig2 . re ’ , STATUS= ’ UNKNOWN ’ ) 
OPEN ( 3 , F I LE= ’ Ine ig 1 . im ’ , STATUS= ’ UNKNOWN ’ ) 
OPEN (4 , FILE= ’ lneig2 . im ’ , STATUS= ’UNKNOWN ’ ) 
OPEN (7 , FILE= ’ IMSLl . re ’ , STATUS= ’ UNKNOWN ’ ) 
OPEN ( 8 , FILE= ’ IMSL2 . re ’ , STATUS= ’ UNKNOWN ’ ) 
OPEN (9 , FILE= ’ IMSL 1 . im ’ , STATUS= ’UNKNOWN ’ ) 
OPEN ( 10 , FILE= ’ IMSL2 . im ’ , STATUS= ’ UNKNOWN ’ ) 
OPEN (20 , FILE= ’ disci ’ , STATUS= ’ UNKNOWN ’ ) 
OPEN (23 , FILE= ’ phase ’ , STATUS= ’ UNKNOWN ’ ) 

C  END  CHECKOUT  CODE 


c 
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c 

c  initialize  at  t=0,  read  if  t  not  zero . 

c 

if(t  .eq.  O.dO)  then 
c 

c  beginning  a  run . read  state 

c 

read  (*,*)  x(l , 1) ,x(2 , 1) 
c 

c  initialize  Phi  =  I,  integrals  of  eigenvlaues  to  zero 

c 

do  100  i  =  3,10 

x(i,l)  =  O.dO 

100  continue 

do  101  i  =  3,6,3 

x(i,l)  =  l.dO 

101  continue 

do  102  i  =  1,2 

eig(i,l)  =  (O.dO,  O.dO) 
npi(i)  =  0 

102  continue 

c  initialize  swap  buffer 

do  103  i  =  1,2 

iswap(i)  =  i 

103  continue 

c  mode  1,  no  bif  checking  at  first 

mode  =  1 

else 

c 

c  resuming  a  run. . . . 

c 

do  104  i  =  1,9,2 

read  (*,*)  x(i , 1) ,x(i+l ,  1) 

104  continue 

c  read  discriminants 

read  (*,*)  discl(4) 
c  read  last  eigenvalues 

do  105  i  =  1,4 

read  (*,*)  wpp(l) ,wpp(2) 

eig(i,4)  =  wp 

105  continue 

read  (*,♦)  imode 
read  (*,*)  (npi(j j) , j j=l ,2) 
read  (*,*)  (iswapCj j ) , j j=l ,2) 
c  begin  in  normal  mode 
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mode  =  0 

endlf 

read  (*,*)  tmax,nstp,nskp 
read  (*,*)  idebug 
c 

nn  =  10 

hh  =  (tmax-t)/(dble(nstp*nskp)) 
nxt  =  0 

call  haming(nxt) 

C  DEBUG 

c  write  (13,*)  ’eig  before  phase  call’ 
c  do  254  i  =  1,2 
c  write  (13,*)  ’eig’,i 

c  do  254  j  =  1,2 

c  write  (13,*)  j,eig(i,j) 

c  254  continue 
C  END 

if (nxt  .eq.  0)  stop  77 

c  adjust  startup  phases 

do  29  i  =  2,4 

call  phased) 

29  continue 
C  DEBUG 

c  write  (13,*)  ’eig  after  phase  calls’ 
c  do  255  i  =  1,2 

c  write  (13,*)  ’eig’,i 

c  do  255  j  =  1,2 

c  write  (13,*)  j,eig(i,j) 

c  255  continue 
C  END 

write  (13,*)  ’haming  initialized’ 
c  cycle  5  times  before  bifurcation  checking 

do  30  i  =  1,5 

call  haming(nxt) 
call  phase(nxt) 

30  continue 

c  turn  on  bifurcation  checking 

mode  =  0 
c 

c  open  output  files 

c 
c 
c 

c  integration  loop 

c 
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do  1000  istp  =  l.nstp 
c 

c  do  eigenvalue  calc  every  other  time.... 

c 

if (mod(istp,l)  .eq.  0)  then 
do  500  i  =  1,2 

do  500  j  =  1,2 

phi(i,j)  =  x(2*j+i,nxt) 

500  continue 

c 

c  call  eigrf (phi, 2,2,1, w,vec, 2, work, ier) 

c 

c  the  eigrf  subroutine  was  not  in  the  IMSL  library  for 

c  eigenvalue  computation,  but  devcrg  is 

c 

call  devcrg(2,phi,2,w,vec,2) 
c 

do  600  i  =  1,2 
wp  =  w(i) 

c  calc  complex  log  w(i) 

c 

if(t  .eq.  O.dO)  go  to  600 

omegr  =  dlog(dsqrt(  wpp(l)*wpp(l)  + 

1  wpp(2)*wpp(2)  )) 

omegi  =  datan2(  wpp(2)  ,  wpp(l)  ) 
c 

c  the  following  if  statement  is  to  write  the  log  mag  and  phase  of  the 

c  IMSL  generated  eigenvalues  in  the  same  order  as  our  literal  method 

c  (i.e.  the  calceig  output):  largest  to  smallest,  or  positive  complex 

c  to  negative  complex 
c 

if  (i.eq.  Dthen 
write  (8,9)  t, omegr 
write  (10,9)  t, omegi 
elseif ( i . eq . 2) then 
write  (7,9)  t, omegr 
write  (9,9)  t, omegi 
endif 

600  continue 

endif 

9  format(2x,2(e20. 13,lx)) 

C  CHECKOUT  CODE:  write  the  "literal"  method  log  mag  and  phase  into 
file  eigKor  2).re(im) 
do  700  i  =  1,2 

wp  =  eig(i,nxt) 
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c 


write  (i,9)  t,wpp(l) 
write  (i+2,9)  t,wpp(2) 
700  continue 

C  END  CHECKOUT  CODE 


c  write  eigenvalues  from  phi  by  "literal"  method 

c  over  t,  as  a  function  of  inverse  time 

c 

C  PRODUCTION  CODE 

if(t  .ne.  O.dO)  then 

8  format(4(lx,el8. 10) ,/ ,lx,el8. 10) 

write  (11,8)  l.d0/t,x(7,nxt)/t,x(8,nxt)/t 
write  (12,8)  1 .dO/t ,x(9,nxt)/t ,x(10,nxt)/t 
endif 

C  END  PRODUCTION  CODE 
c 

c  do  the  integration 

c 

do  900  jskp  =  l,nskp 
c 

call  haming(nxt) 
c 

call  split (nxt) 
c 

c  adjust  phase  after  bif  checks 

call  phase (nxt) 
c 

900  continue 

c 

1000  continue 
c 

c  dump  final  state 
c 

call  dump (nxt) 
c 

stop 

and 


c$INCLUDE 
c$INCLUDE 
c$ INCLUDE 
clINCLUDE 
c$ INCLUDE 
c$INCLUDE 
c$ INCLUDE 


’haming.for’ 
’rhs_vpl .for’ 
’calceig.vpl .for’ 
’split_vpl .for’ 
’split2_vpl .for’ 
’checks.vpl .for’ 

’  dump. vpl. for’ 
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c$INCLUDE:  ’phase.vpl .for' 
c$INCLUDE:  ’ amat.vpl .for ' 
c 

include  'haming.vpl .for’ 
include  'rhs.vpl .for’ 
include  ’calceig_vpl.for’ 
include  ’split.vpl .for’ 
include  ’split2_vpl .for’ 
include  ’checks.vpl .for’ 
include  ’ dump _vpl .for’ 
include  ’phase_vpl .for’ 
include  ’ amat.vpl .for ’ 
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c 

subroutine  haming(nxt) 
c 

c  haming  is  an  ordinary  differential  equations  integrator 
c  it  is  a  fourth  order  predictor-corrector  algorithm 
c  which  means  that  it  carries  along  the  last  four 
c  values  of  the  state  vector,  and  extrapolates  these 

c  values  to  obtain  the  next  value  (the  prediction  part) 

c  and  then  corrects  the  extrapolated  value  to  find  a 

c  new  value  for  the  state  vector, 

c 

c  the  value  nxt  in  the  call  specifies  which  of  the  4  values 
c  of  the  state  vector  is  the  "next"  one. 

c  nxt  is  updated  by  haming  automatically,  and  is  zero  on 

c  the  first  call 

c 

c  the  user  supplies  an  external  routine  rhs(nxt)  which 
c  evaluates  the  equations  of  motion 

c 

common  /ham/  x,y(10,4) ,f (10,4) ,errest(10) ,n,h,mode,e 
double  precision  x,y,f ,errest,h,hh,xo,tol,e 
c 

c  all  of  the  good  stuff  is  in  this  common  block, 

c  X  is  the  independent  variable  (  time  ) 

c  y(6,4)  is  the  state  vector-  4  copies  of  it,  with  nxt 

c  pointing  at  the  next  one 

c  f(6,4)  are  the  equations  of  motion,  again  four  copies 
c  a  call  to  rhs(nxt)  updates  an  entry  in  f 

c  errest  is  an  estimate  of  the  truncation  error  -  normally  not 
c  used 

c  n  is  the  number  of  equations  being  integrated  -  6  or  42  here 

c  h  is  the  time  step 

c  mode  is  0  for  just  EOM,  1  for  both  EOM  and  EOV 

c 

tol  =  O.OOOOOOOOOld+00 

c  switch  on  starting  algorithm  or  normal  propagation 

if (nxt)  190,10,200 
c 

c  this  is  hamings  starting  algorithm. ...  a  predictor  -  corrector 

c  needs  4  values  of  the  state  vector,  and  you  only  have  one-  the 

c  initial  conditions. 

c  haming  uses  a  Picard  iteration  (slow  and  painfull)  to  get  the 

c  other  three. 

c  if  it  fails,  nxt  will  still  be  zero  upon  exit,  otherwise 

c  nxt  will  be  1,  zoid  you  are  all  set  to  go 


A-8 


c 


10  xo  =  X 

hh  =  h/2.0d+00 
call  rhs(l) 
do  40  1  =  2,4 

X  =  X  +  hh 

do  20  i  =  l,n 

20  y(i,l)  =  yCi,l-l)  +  hh*f(i.l-l) 
call  rhs(l) 

X  =  X  +  hh 

do  30  i  =  l,n 

30  y(i,l)  =  y(i,l-l)  +  h*f(i,l) 

40  call  rhs(l) 
jsw  =  -20 
50  isw  =  1 

do  120  i  =  l,n 

hh  =  y(i,l)  +  h*(  9.0d+00*f (i,l)  +  19.0d+00*f (i,2) 

1  -  5.0d+00*f (i,3)  +  f(i,4)  )  /  24.0d+00 

if(  dabs(  hh  -  y(i,2))  .It.  tol  )  go  to  70 
isw  =  0 

70  y(i.2)  =  hh 

hh  =  y(i,l)  +  h*(  f(i,l)  +  4.0d+00*f (i,2)  +  f (i,3))/3.0d+00 
if(  dabs(  hh-y(i,3))  .It.  tol  )  go  to  90 
isw  =  0 

90  y(i,3)  =  hh 

hh  =  y(i,l)  +  h*(  3.0d+00*f (i,l)  +  9.0d+00*f (i,2)  +  9.0d+00*f (i,3) 
1  +  3.0d+00*f (i,4)  )  /  8.0d+00 

if(  dabs(hh-y(i,4))  .It.  tol  )  go  to  110 
isw  =  0 

110  y(i,4)  =  hh 
120  continue 
X  =  xo 

do  130  1  =  2.4 

X  =  X  +  h 
130  call  rhs(l) 

if(isw)  140,140,150 
140  jsw  =  jsw  +  1 

if(jsw)  50,280,280 
150  X  =  xo 
isw  =  1 
jsw  =  1 

do  160  i  =  l,n 
160  errest(i)  =  0.0 
nxt  =  1 
go  to  280 
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190  jsw  =  2 

nxt  =  iabs(nxt) 
c 

c  this  is  hamings  normal  propagation  loop  - 
c 

200  X  =  X  +  h 

npl  =  mod(nxt,4)  +  1 
go  to  (210,230) ,isw 
c  permute  the  index  nxt  modulo  4 
210  go  to  (270,270,270,220) , nxt 
220  isw  =  2 

230  nm2  =  mod(npl,4)  +  1 

nml  =  mod(nm2,4)  +  1 

npo  =  mod(nml,4)  +  1 

c 

c  this  is  the  predictor  part 
c 

do  240  i  =  l,n 

f(i,nm2)  =  y(i,npl)  +  4.0d+00*h*(  2.0d+00*f (i,npo)  -  f(i,nml) 

1  +  2.0d+00*f (i,nm2)  )  /  3.0d+00 

240  y(i,npl)  =  f(i,nm2)  -  0.925619835d0*errest(i) 
c 

c  now  the  corrector  -  fix  up  the  extrapolated  state 
c  based  on  the  better  value  of  the  equations  of  motion 
c 

call  rhs(npl) 
do  250  i  =  l,n 

y(i,npl)  =  (  9 .0d+00*y(i ,npo)  -  y(i,nm2)  +  3.0d+00*h*(  f(i,npl) 
1  +  2.0d+00*f (i,npo)  -  f(i,nml)  )  )  /  8.0d+00 

errest(i)  =  f(i,nm2)  -  y(i,npl) 

250  y(i,npl)  =  y(i,npl)  +  0 .074380 1653d0  *  errest(i) 
go  to  (260,270) , jsw 
260  call  rhs(npl) 

270  nxt  =  npl 
280  return 
end 
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c 

c 

subroutine  rhs(nxt) 
c 

c  Van  der  Pol’s  "simple"  equation 

c  eom  and  phi  matrix,  integrals  of  eigenvalues  dt 

c 

common  /ham/t,x(10,4) ,f (10,4) ,err(10) ,nn,hh,mode ,e 
double  precision  t,x,f ,err,hh,e 

common  /eval/  eig(2,4) ,npi(2) ,discl(4) ,dtbif ,iswap(2) 
complex* 16  eig.wp 

double  precision  disci ,dtbif ,wpp(2) 
common  /debug/  idebug 
double  precision  a(2,2) 
double  precision  xx(2) 
equivalence  (wp,wpp) 
c 

C  DEBUG 

c  write  (*,*)  ’rhs,  nxt=’,nxt 

c  write  (*,*)  ’timestep’,  hh 

c  write  (*,*)  ’full  state/rate  dump  at  entry’ 
c  do  333  i  =  l,nn 

c  write  (*,*)  i,x(i,nxt) ,f (i,nxt) 

c  333  continue 
C  END 
c 

c  extract  state 

c 

do  1  i  =  1 , 2 

xx(i)  =  x(i,nxt) 

1  continue 
c 

c  Van  der  Pohl’s  eom 

c 

f(l,nxt)  =  x(2,nxt) 

f (2,nxt)  =-x(l ,nxt)+(e-e*x(l ,nxt)*x(l ,nxt))*x(2,nxt) 
c 

c  A(t)  matrix 

c 

call  amat(nxt,a) 
c 

c  phi  dot  =  a  phi 

c  phi  stored  by  cols 

c 

do  50  i  =  1,2 


All 


do  50  j  =  1,2 

f(i+2*j,nxt)  =  O.dO 
do  50  k  =  1,2 

f(i+2*j,nxt)  =  f(i+2*j,nxt)  +  a(i,k)*x(k+2*j ,nxt) 

50  continue 
c 

c  calculate  eigenvalues 

c 

if(  t  .ne.  0)  then 

call  caleig(nxt) 

endif 

c 

c  eigenvalue  average  eom 

c 

if(  t  .ne.  0)  then 
wp  =  eig(l,nxt) 
f(7,nxt)  =  wpp(l)/t 
f(9,nxt)  =  wpp(2)/t 
wp  =  eig(2,nxt) 
f(8,nxt)  =  wpp(l)/t 
f(10,nxt)  =  wpp(2)/t 

else 

f(7,nxt)  =  O.dO 
f(8,nxt)  =  O.dO 
f(9,nxt)  =  O.dO 
f(10,nxt)  =  O.dO 

endif 

c 

C  DEBUG 

c  ifCidebug  .ge.  5.)  then 

c  write  (13,*)  ’rhs,  full  state/rate  dump,  nxt=’,nxt 

c  do  666  i  =  1,10 

c  write  (13,*)  x(i ,nxt) ,f (i ,nxt) 

c  666  continue 

c  endif 

C  END 
c 

return 

end 


c 

c 

subroutine  amat(nxt,a) 
c 

c  calculate  A  matrix  for  Van  der  Pol’s  ’simple  ’  system, 
c 

common  /ham/  t,x(10,4) ,f (10,4) ,err(10) ,nn,hh,mode,e 
double  precision  t,x,f ,err,hh,e 
double  precision  xx(2),a(2,2) 
c 
c 

c  extract  state 
c 

do  5  i  =  1 , 2 

xx(i)  =  x(i,nxt) 

5  continue 
c 
c 

a(l,l)=0.d0 

a(l,2)=l.d0 

a(2, 1)=-1 .d0-2.d0*e*xx(l)*xx(2) 
a(2,2)=e-e*xx(l)*xx(l) 
c 

return 

end 
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c 

c 

subroutine  caleig(nxt) 
c 

c  calculate  eigenvalues  of  Phi  matrix  for  simple  2D  van  der  Pol 
c  system,  from  integrated  Phi  matrix  elements 
c 

common  /ham/  t,x(10,4) ,f (10,4) ,err(10) ,nn,hh,mode,e 
double  precision  t,x,f ,err,hh,e 

common  /eval/  eig(2,4) ,npi(2) ,discl(4) ,dtbif ,iswap(2) 
complex* 16  eig 

double  precision  disci ,dtbif ,a,b,c 
complex*16  fig(4),ca,cb 
common  /debug/  idebug 
c 

double  precision  phi(2,2) ,wpp(2) ,wqq(2) ,vmag 
double  precision  twopi 
complex*16  wp,wq,cdisc 
equivalence  (wp ,wpp) , (wq,wqq) 

C  DEBUG 

c  if(idebug  .ge.  2)  write  (13,*)  ’enter  calceig’ 

C  END 
c 

twopi  =  2,d0*3.141592653589793d0 
nml  =  nxt  +  3 

if(nml  .gt.  4)  nml  =  nml  -  4 
c 

c  extract  phi,  stored  by  cols 
c 

do  10  i  =  1,2 

do  10  j  =  1,2 

phi(i,j)  =  x(2*j+i,nxt) 

10  continue 
c 

c  solve  for  eigenvalues  of  phi  (2x2),  and  discriminant  (disci) 

c  disci  is  the  argument  under  the  square  root  in  the  quadratic  eq 

c  solution.  Therefore,  it  indicates  whether  we  have  a  real,  or 

c  complex  eigenvalue, 

c 

b  =  -phi(l,l)-phi(2,2) 

c  =  (phi(l,l)*phi(2,2))  -  (phi(i ,2)*phi(2 , 1) ) 
c 

C  DEBUG 

c  if (idebug  .ge.  2)  write  (13,*)  ’A’,b 

C  END 
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c 

disci (nxt)  =  .25d0*b*b  -  c 
c 

c  note:  tl-e  following  write  statement  to  file  20,  called  disc,  is 

c  used  to  plot  the  discriminant  vs  time.  This  enables  us  to  see 

c  the  long  term  behavior  of  the  discriminant,  and  therefore  correct 

c  for  bifurcations 

c 

write(20,20)  t,discl(nxt) 

20  format (2x ,2(e20. 13 , lx) ) 

c 
c 
c 

if  (discl(nxt)  .ge.  O.dO)  then 
wpp(l)  =  dsqrt  (discl(nxt)) 
wpp(2)  =  O.dO 

else 

wpp ( 1 )  =  O.dO 

wpp(2)  =  dsqrt  (dabs  (disci (nxt)  )  ) 
endif 
c 

wqq(l)  =  -b/2.d0 

wqq(2)  =  O.dO 

c 

C  DEBUG 

c  if(idebug  .ge.  2)  then 

c  write  (13,*)  ’disc’ ,discl(nxt) 

c  write  (13,*)  ’sqrt’,wp 

c  endif 

C  ENDIF 
c 

fig(l)  =  wq  +  wp 
fig(2)  =  wq  -  wp 
c 

c  iswap  will  keep  track  of  bifurcation  "direction" 

c 

do  200  i  =  1,2 

eig(i,nxt)  =  f ig(iswap(i) ) 

200  continue 

c 

c  calculate  logs  of  eigenvalues 

c 

do  100  i  =  1,2 

wp  =  eig(i,nxt) 

wqq(l)  =  dlog(  dsqrt(  vpp(l)*wpp(l)  +  wpp(2)*wpp(2)  )) 
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wqq(2)  =  datan2(  wpp(2),  wpp(l)  )  +  dble(npi(i))*twopi 

C  DEBUG 

c  if(idebug  .ge.  2)  then 

c  write  (13,*)  ’t=’,t 

c  write  (13,*)  ’i,  lambda’ ,i,eig(i ,nxt) 

c  write  (13,*)  ’i.  In  lam’,i,wq 

c  endif 

C  END 

eig(i,nxt)  =  wq 
100  continue 
c 

C  DEBUG 

c  if(idebug  .ge.  2)  write  (13,*)  ’exit  calceig’ 

C  END 

return 

end 
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c 

c 

subroutine  split (nxt) 
c 

c  check  for  root  bifurcations  and  fix  up  results 
c 

common  /ham/  t,x(10,4) ,f (10,4) ,err(10) ,nn,hh,mode,e 
double  precision  t,x,f ,err,hh,e 

common  /eval/  eig(2,4) ,npi(2) ,discl(4) ,dtbif ,iswap(2) 
complex* 16  eig.etempl ,etemp2 
double  precision  disci ,dtbif .dtempl ,dtemp2 
common  /debug/  idebug 
c 

double  precision  wpp(2) ,wqq(2) 
complex* 16  wp.wq 
equivalence  (wp.wpp) , (wq.wqq) 
c 

c  do  bifurcation  checks?  not  on  startup 

c 

if (mode  .ne.  0)  return 
c 
c 

c  BIFURCATION  CHECKS! 

c 

nml  =  nxt  +  3 

if(nml  .gt.  4)  nml  =  nml  -  4 
c 

c  check  for  isolated  bifurcation 

c 

if(  discl(nxt)*discl(nml)  .It.  O.dO  )  then 
c 

c  DEBUG 
c 
c 
c 
c 

c  END 

c 
c 
c 

C  DEBUG 


if (idebug  .gt.  0)  then 

write  (13,*)  ’isolated  bif’ 

write  (13,*)  ’real  pts ’ ,wpp(l) ,wqq(l) 

endif 


call  split2(nxt) 

do  reasonableness  checks,  write  breakpt  file 

temp  higher  debug  level .... 
ida  =  idebug 

if (idebug  .ge.  1)  idebug  =  2 
if (ida  .gt.  0)  then 


A- 17 


C  END 


write  (13,*)  ’post  split2  checks’ 
endif 

call  checks (nxt) 

C  DEBUG  reset  debug  level 

idebug  =  ida 

C  END 
c 

c  bifurcation  introduces  a  discontinuity  in  integral  slopes 

c  haming  must  be  restarted 

c 

do  200  i  =  1,10 

x(i,l)  =  x(i,nxt) 

200  continue 

etempl  =  eig(l,nml) 
etemp2  =  eig(l,nxt) 
eig(l,l)  *  etemp2 
eig(l,4)  =5  etempl 
etempl  =  eig(2,nml) 
etemp2  =  eig(2,nxt) 
eig(2,l)  =  etemp2 
eig(2,4)  =  etempl 
dtempl  =  disci (nml) 
dterap2  *  disci (nxt) 
disci (1)  *  dtemp2 
disci (4)  ®  dtempl 
nxt  =  0 

C  DEBUG 

c  ida  =  idebug 

c  idebug  =  6 

C  END 

call  haming(nxt) 
if (nxt  ,eq.  0)  then 

write  (13,*)  ’haming  restart  failure  after  bifurcation’ 
stop 
endif 

C  DEBUG 

c  idebug  =  ida 

C  END 

c  PRAY  HERE  THAT  ANOTHER  BIF  DIDN’T  OCCUR 

c 

endif 

c 

C  DEBUG 

c  write  (13,*)  ’exit  split’ 
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C  END 


return 

end 
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c 

c 

subroutine  split2(nxt) 
c 

c  hauidles  isolated  bifurcations  of  a  single  +-  root  pa 
c 

common  /ham/  t,x(10,4) ,f (10,4) ,err(10) ,nn,hh,mode,e 
double  precision  t,x,f ,err,hh,e 

common  /eval/  eig(2,4) ,npi(2) ,discl(4) .dtbif ,iswap(2) 
complex* 16  eig 

double  precision  disci, dtbif 
common  /debug/  idebug 
complex* 16  wp,wq 

double  precision  dnxt ,dnml ,dnm2,dnm3 
double  precision  wpp(2) ,wqq(2) ,ttfopi 
double  precision  p,dp,capd,capdot 

double  precision  sml,sm2,sm3,cml,cm2,cm3,ints(4) ,splus2,tbif 
complex* 16  eig0(2) 
equivalence  (wp ,wpp) , (wq,wqq) 
c 

twopi  =  2.d0*3.141592653589793d0 


C  DEBUG 

if (idebug  .ne.  0)  then 
write  (13,*)  ’  ’ 

write  (13,*)  ’isolated  bifurcation,  root  pair’ 
write  (13,*)  ’near  t=’,t 


END 


endif 

set  indices 

to 

extract  roots 

nml  =  nxt  + 

3 

if (nml 

.gt. 

4) 

nml 

=  nml 

-  4 

nm2  =  nxt  + 

2 

if (nm2 

•gt. 

4) 

nm2 

=  nm2 

-  4 

nm3  =  nxt  + 

1 

if  (nm3 

•  gt. 

4) 

nm3 

=  nm3 

-  4 

il  =  1 

12  =  2 

dnxt  = 

disci (nxt) 

dnml  = 

discl(nml) 

dnm2  = 

disci (nm2) 

dnm3  = 

disci (nm3) 
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C  DEBUG 

if(idebug  .ne.  0)  then 

write  (13,*)  ’roots  ’,il,i2 
write  (13,*)  ’dnxt’ ,dnxt 
write  (13,*)  ’dnml ’ ,dnml 
write  (13,*)  ’dnm2’,dnm2 
write  (13,*)  ’ dnm3 ’ , dnm3 
endif 

C  END 
c 

c  determine  bifurcation  time 
c 

p  =  -dnxt  /(  dnxt  -  dnml  ) 
do  50  i  =  1,15 

capd  =  (p+1 .d0)*(p+2.d0)*(p+3.d0)*dnxt/6.d0 

1  -  p*(p+2.d0)*(p+3.d0)*dnml/2.d0 

2  +  p*(p+l .d0)*(p+3.d0)*dnm2/2.d0 

3  -  p*(p+l .d0)*(p+2.d0)*dnm3/6.d0 

capdot  =  ((p+2.d0)*(p+3.d0)+(p+l.d0)*(p+3.d0)+ 

1  (p+1 .d0)*(p+2.d0))*dnxt/6.d0 

2  -  ((p+2.dO)*(p+3.dO)+p*(p+3.dO)+p*(p+2.dO)) 

3  *  dnml/2.d0 

4  +  ((p+1 .d0)*(p+3.d0)+p*(p+3.d0)+p*(p+l .dO)) 

5  *  dnm2/2.d0 

6  -  ((p+l.d0)*(p+2.d0)+p*(p+2.d0)+p*(p+l.d0)) 

7  *  dnm3/6.d0 

dp  =  -capd/capdot 

c  DEBUG 

ifddebug  .ne.  0)  then 

write  (13,*)  ’D,  Ddot,  p,  dp’ , capd, capdot, p, dp 
endif 

c  END 

p  =  p  +  dp 

if (dabs (dp)  .It.  l.d-7)  go  to  55 
50  continue 

write  (13,*)  ’SPLIT2;  bifurcation  time  iteration  failed’ 
stop  63 
55  continue 
dtbif  =  p*hh 
C  DEBUG 

if(idebug  .ne.  0)  then 

write  (13,*)  ’dtbif ’, dtbif 
endif 

C  END 
c 
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c  determine  bifurcation  type 

c 

if(dnxt  .It.  O.dO)  go  to  2000 
c 

1000  continue 
c 

C  *  ♦ 

c  *  current  discriminant  positive,  imag  pair  ->  reals  * 

c  *  * 

C 

C  DEBUG 

if(idebug  .ne.  0)  then 

write  (13,*)  ’imag  ->  real’ 
end  if 

C  END 
c 

itype  =  1 
c 

c  reset  two  pi  multipliers  and  correct  new  roots 
c 

call  phase(nxt) 
c 

c  now  fix  up  average  eigenvalue  integrals  over  bifurcation 
c 

go  to  5000 
c 

2000  continue 
c 

c  ♦♦♦♦♦********♦*♦**♦**♦*♦♦♦*♦********♦*♦*♦♦*♦***♦***♦♦**♦♦♦ 

c  *  * 

c  *  current  discriminant  neg,  real  pair  ->  imaginaries  * 

c  *  * 

c  ♦♦*♦♦♦♦♦****♦♦**♦♦****♦*♦*♦*****♦♦*♦♦♦♦♦*>(■♦♦♦♦*♦♦*+**♦♦♦** 

c 

C  DEBUG 

ifCidebug  .ne.  0)  then 

write  (13,*)  ’real  ->  imag’ 
endif 

C  END 
c 

itype  =  2 
c 

c  swap  this  pair  as  necessary  to  insure  that  root  with  positive 
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c  imag  part  always  increases . 

c 

c  who  has  positive  imag  paurt?  check  last  root 
ipos  =  il 
wp  =  eig(il,nml) 
if(wpp(2)  .It.  O.dO)  ipos  =  i2 

c  bifurcation  on  left  or  right  side  of  zero?  Check  root  phase 

c  prior  root  is  real,  and  prior  phase  is  a  multiple  of  pi... 

wp  =  eig(ipos,nml) 
c  reduce  mod  pi 

icyc  =  idnintC  2.d0*wpp(2)/twopi  ) 
iside  =  +1 

if(  modCicyc,  2)  .ne.  0)  iside  =  -1 

c  reduce  NEW  root  phase  mod  2  pi  to  see  if  increasing  or  decreasing 
wp  =  eigCipos ,nxt) 
jcyc  =  idint(  wpp(2)/twopi  ) 
wpp(2)  =  wpp(2)  -  dble(jcyc)*twopi 
C  DEBUG 

if(idebug  .gt.  0)  then 

write  (13,*)  ’positive  root  now’ , ipos, eig(ipos ,nxt) 
write  (13,*)  ’positive  root  pre’ , ipos ,eig( ipos ,nml) 
write  (13,*)  ’prior  phase  /  pi’, icyc 
write  (13,*)  ’new  reduced  phase* ,wpp(2) 
write  (13,*)  ’axis  side’, iside 

endif 

C  END 

c  determine  if  swapping  needed. . .side  dependent 
if (iside  .gt.  0)  then 

if(wpp(2)  .gt.  4.71d0)  go  to  2005 

else 

if(wpp(2)  .It.  3.141592653589d0)  go  to  2005 

endif 

c  no  swapping  needed 

C  DEBUG 

if(idebug  .ne.  0)  write  (13,*)  ’no  swapping  needed’ 

C  END 

go  to  2010 
2005  continue 
c 

c  swap  this  root  pair! 

c 

ineg  =  il 

if(ineg  .eq.  ipos)  ineg  =  i2 
wp  =  eig(ipos ,nxt) 

wpp(2)  =  wpp(2)  +  (dble(npi(ineg))  -  dble(npi(ipos)))*twopi 
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wq  =  eig(ineg,nxt) 

wqq(2)  =  wqq(2)  +  (dble(npi(ipos) )  -  dble(npi(ineg)))*twopi 
eigCipos ,nxt)  =  wq 
eig(ineg,nxt)  =  wp 

c  swap  iswap  order.... its  tied  to  calceig  order 
isw  =  iswap (ipos) 
iswap(ipos)  =  iswap(ineg) 
iswap (ineg)  =  isw 
C  DEBUG 

if(idebug  .ne.  0)  then 

write  (13,*)  ’roots  swapped’ 

endif 

C  END 
c 

2010  continue 
c 

c  reset  two  pi  multipliers  and  correct  new  roots 

c 

call  phase (nxt) 
c 

C  4<4:*********************************<|:*#**:ti4:*^!|:><<!|i!4c«4c:|i:|c4i 

c  *  carry  average  eigenvalue  ints  over  bif  * 

C  :iHilttt*********************************!**************^* 

c 

5000  continue 
C  DEBUG 

if(idebug  .ne.  0)  then 

write  (13,*)  ’pre-bif  eigenvalues’ 

write  (13,*)  eig(l,nml) 

write  (13,*)  eig(2,nml) 

write  (13,*)  ’post-bif  eigenvalues’ 

write  (13,*)  eig(l,nxt) 

write  (13,*)  eig(2,nxt) 

endif 

C  END 
c 

c  extrapolate  integrals,  log  eigenvalues  into  bif,  s  interp 

c 

c  s  factors,  (dtbif  is  negative) 

sml  =  dsqrt(  dabs(  hh  +  dtbif  )) 

sm2  =  dsqrt(  dabs(  hh  +  dtbif  )  +  hh  ) 

sm3  =  dsqrt(  dabs(  hh  +  dtbif  )  +  2.d0*hh  ) 

c 

c  extrap  coef 

cml  =  sm2*sm3/(  (sml-sm2)*(sml-sm3)  ) 
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cm2  =  sml*sm3/(  (sm2-sml)*(sm2-sm3)  ) 
cm3  =  sml*sm2/(  (sm3-sml)*(sm3-sm2)  ) 
c 

GigO(l)  =  cml*eig(l ,nml)  +  cm2*eig(l ,nm2)  +  cm3*eig(l ,nm3) 
eig0(2)  =  cml*eig(2,nml)  +  cm2*eig(2,nm2)  +  cm3*eig(2 ,nm3) 
c  DEBUG 

c  ifCidebug  .ne.  0)  then 

c  write  (13,*)  ’extrapolation  factors’ 

c  write  (13,*)  cml,cm2,cm3 

c  write  (13,*)  ’eigenvalues  extrapolated  to  s=0’ 

c  write  (13,*)  eigO(l) 

c  write  (13,*)  eig0(2) 

c  endif 

C  END 
c 

ints(l)  =  cml*x(7,nml)  +  cm2*x(7,nm2)  +  cm3*x(7,nm3) 

ints(2)  =  cml*x(8,nml)  +  cm2*x(8,nm2)  +  cm3*x(8,nm3) 

ints(3)  =  cml*x(9,nml)  +  cm2*x(9,nm2)  +  cm3*x(9,nm3) 

ints(4)  =  cml*x(10,nml)  +  cm2*x(10,nm2)  +  cm3*x(10,nm3) 

C  DEBUG 

c  ifddebug  .ne.  0)  then 

c  write  (13,*)  ’last  3  pre  bif  integrals’ 

c  write  (13,*)  ’7  ’ ,x(7,nml) ,x(7,nm2) ,x(7,nm3) 

c  write  (13,*)  ’8  ’ ,x(8,nml) ,x(8,nm2) ,x(8,nm3) 

c  write  (13,*)  ’9  ’ ,x(9,nml) ,x(9,nm2) ,x(9,nm3) 

c  write  (13,*)  ’ 10’ ,x(10,nml) ,x(10,nm2) ,x(10,nm3) 

c  write  (13,*)  ’averaged  integrals  extrap  to  s=0’ 

c  write  (13,*)  ints(l) ,ints(3) 

c  write  (13,*)  ints(2) ,ints(4) 

c  endif 

C  END 
c 

c  integrate  out  of  singularity,  s  integral.  In  eigenvalues 

c  assumed  linear  in  s 

c 

splus2  =  dabs(  dtbif  ) 
tbif  =  t  +  dtbif 
c 

wp  =  eigO(l) 
wq  =  eig(l,nxt) 

x(7,nxt)  =  ints(l)  +  2.d0*splus2*(  wpp(l)/6.d0  +  wqq(l)/3.d0)/tbif 
x(9,nxt)  =  ints(3)  +  2.d0*splus2*(  wpp(2)/6.d0  +  wqq(2)/3.d0)/tbif 
wp  =  eig0(2) 
wq  =  eig(2,nxt) 

x(8,nxt)  =  ints(2)  +  2 .d0*splus2*(  wpp(l)/6.d0  +  wqq(l)/3.d0)/tbif 
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x(10,nxt)=  ints(4)  +  2 .d0*splus2*(  wpp(2)/6.d0  +  wqq(2)/3 .dO)/tbif 
C  DEBUG 

c  if(idebug  .ne.  0)  then 

write  (13,*)  ’ints  extrap  to  end  of  timestep’ 
c  write  (13,*)  y(7,nxt) ,x(9,nxt) 

c  write  (13,*)  x(8,nxt) ,x(10,nxt) 

c  endif 

C  END 


return 

end 
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c 

c 

subroutine  phase (nxt) 
c 

c  fix  phase  on  current  eigenvalues,  update  npi 
c 

common  /ham/  t,x(10,4) ,f (10,4) ,err(10) ,nn,hh,mode,e 
double  precision  t,x,f ,err,hh,e 

common  /eval/  eig(2,4) ,npi(2) ,discl(4) ,dtbif ,iswap(2) 

complex* 16  eig 

double  precision  disci, dtbif 

common  /debug/  idebug 

complex* 16  wp,wq 

double  precision  wpp(2) ,wqq(2) ,ttfopi 
equivalence  (wp,wpp) , (wq.wqq) 
c 

twopi  =  2.d0*3.141592653589793d0 
c 
c 

c  set  indices  to  extract  roots 
c 

nml  =  nxt  +  3 

if(nml  .gt.  4)  nml  =  nml  -  4 
c 

c  reset  two  pi  multipliers  and  correct  new  roots 
c 

do  100  il  =  1,2 

wp  =  eig(il,nxt) 
wq  =  eig(il,nml) 

idel  =  nint(  (wqq(2)  -  wpp(2)) /twopi  ) 

C  DEBUG 

if (idebug  .ge.  2)  then 

write  (23,*)  ’npi ’ ,npi(il) 
write  (23,*)  ’eig  ’,il,’  now’,wp 
write  (23,*)  ’eig  ’,il,’  pre’,wq 
write  (23,*)  ’idel’, idel 

endif 

C  END  note;  file  23  is  called  ’phase’ 
c 

if (idel  .ne.  0)  then 

npi(il)  =  npi(il)  +  idel 

wpp(2)  =  wpp(2)  +  dble(idel)*twopi 

eig(il,nxt)  =  wp 

C  DEBUG 

if (idebug  .ge.  2)  then 
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C  END 


write  (23,*)  ’corrected  eig’.wp 
endif 


endif 

100  continue 

return 

end 
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c 

c 

subroutine  checks (nxt) 
c 

c  do  some  reasonableness  checks  for  van  der  Pol’s  "simple  system" 
c 

common  /ham/  t,x(10,4) ,f (10,4) ,err(10) ,nn,hh,mode,e 
double  precision  t,x,f ,err,hh,e 

common  /eval/  eig(2,4) ,npi(2) ,discl(4) .dtbif ,iswap(2) 
complex* 16  eig 
double  precision  disci, dtbif 
common  /debug/  idebug 
c 

double  precision  wpp(2) ,wqq(2) 
complex* 16  wp,wq 
equivalence  (wp,wpp) , (wq,wqq) 
c 

c  check  eigenvalues  (imaginary  parts)  for  pairedness 
c 

do  100  i  =  1,1 

wp  =  eig(2*(i-l)+l ,nxt) 

wq  =  eig(2*(i-l)+2,nxt) 

c  imag  parts  negs 

if (dabs(wpp(2)+wqq(2) )  .gt.  O.ldO)  then 
itrip  =  2 
go  to  666 

end  if 

100  continue 
c 

c  check  averaged  integrals  for  pairedness 

c 

if(  dabs(x(9,nxt)+x(10,nxt))  .gt.  O.OOOSdO)  then 
itrip  =  4 
go  to  666 

endif 

c 

c  we’ve  escaped  the  big  red  button  one  more  time 

c  save  here . 

c 

write  (13,*)  ’breakpoint  file  written  at  t=’,t 
call  dump (nxt) 


c  ABORT  PROCESSING 
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c 


666  continue 
c 

write  (13,*)  ’CHECK  ABORT’ 
write  (13,*)  ’at  t=’ ,t 
c 

if(itrip  .eq.  1)  then 

write  (13,*)  ’real  parts  not  negs’ 
elseif(itrip  .eq.  2)  then 

write  (13,*)  ’imag  parts  not  negs’ 
elseifdtrip  .eq.  3)  then 

write  (13,*)  ’real  parts  integrals  not  neg’ 

else 

write  (13,*)  ’imag  parts  integrals  not  neg’ 

endif 

c 

nml  =  nxt  -  1 

if(nml  .le.  0)  nml  =  nml  +  4 
nm2  =  nxt  -  2 

if(nm2  .le.  0)  nm2  =  nm2  +  4 
nm3  =  nxt  -  3 

if(nm3  .le.  0)  nm3  =  nm3  +  4 
c 

do  10  j  =  1,2 

write  (13,*)  ’eig’,j 
write  (13,*)  ’  nxt’ ,eig(j ,nxt) 

write  (13,*)  ’  nml ’ ,eig(j ,nml) 

write  (13,*)  ’  nm2’ ,eig(j ,nm2) 

write  (13,*)  ’  nm3’ ,eig(j ,nm3) 

10  continue 
c 
c 

write  (13,*)  ’disci’ 
write  (13,*)  ’  nxt’ , disci (nxt) 

write  (13,*)  ’  nml ’, disci (nml) 

write  (13,*)  ’  iun2’ ,discl(nm2) 

write  (13,*)  ’  nm3’ , disci (nm3) 

c 

stop 

end 
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subroutine  dump(nxt) 


dump  current  state  to  breakpoint  file 

common  /ham/  t,x(10,4) ,f (10,4) ,err(10) ,nn,hh,mode,e 
double  precision  t,x,f ,err,hh,e 

common  /eval/  eig(2,4) ,npi(2) ,discl(4) .dtbif ,iswap(2) 

complex* 16  eig 

double  precision  disci, dtbif 

common  /debug/  idebug 

double  precision  wpp(2) 

complex* 16  wp 

equivalence (  wp,wpp) 

dump  current  state  to  file 

open(15,FILE=’break.pt ’ ) 

123  format(lx,e20.13,lx,e20.13) 
write  (15,123)  t 

do  124  i  =  1,9,2 

write  (15,123)  x(i,nxt) ,x(i+l,nxt) 

124  continue 

125  format(lx,4(i5,lx)) 
write  (15,123)  disci (nxt) 
do  127  i  =  1,2 

wp  =  eig(i,nxt) 

write  (15,123)  wpp(l) ,wpp(2) 

127  continue 

write  (15,125)  (npi(j j) , jj=l,2) 
write  (15,125)  (iswap(j j) , j j=l,2) 

close(15) 


return 

end 


Appendix  B.  Computer  Programs  Used  for  Fast  Fourier 
Transform/ Power  Spectral  Density  Analysis 
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c 

program  fftmain 
c 

c  2D  system 

c  integrate  the  eom  and  variational  equations  for  a  2D  system 
c 

common  /ham/  t ,x(10 ,4) ,f (10 ,4) ,err(10) ,nn,hh,mode,e 
double  precision  t,x,f ,err,hh,e 
common  / debug/  idebug 
double  precision  tmax,dmag 
double  precision  a(2,2) ,work(40) 
c 

c  input  parameters,  Ic’s 

c 

read  (*,*)  t 
read  (*,*)  e 
c 

c  open  output  file 

c 

open ( 13 , FILE= ’ run . out ’ , STATUS= 'UNKNOWN ’ ) 

C  PRODUCTION  CODE 

OPEN ( 1 , FILE= ’ stat e ’ , STATUS= ’ UNKNOWN ’ ) 

OPEN (2 .FILE* ’ varl ’ .STATUS* 'UNKNOWN' ) 

OPEN (3 , FILE* ' var2 ' , STATUS* 'UNKNOWN ' ) 

C  END  PRODUCTION  CODE 

c 

c 

c  initialize  at  t=0, 

c 

c  beginning  a  run . read  state 

c 

read  (*,*)  x(l , 1) ,x(2,l) ,x(3,l) ,x(4, 1) 
c 

c  read  integration  parameters:  max  time,  #  of  steps  in  max  time, 
c  how  often  integration  is  reported 

read  (*,*)  tmax,nstp,nskp 
c 

c  nn  is  the  niunber  of  equations  to  be  integrated,  hh  will  determine 

c  the  time  step  that  the  integrator,  haming,  will  use.  Four  copies 

c  of  the  state  vector  are  kept  around,  thus  the  incrementor  nxt. 
c 

nn  =  4 

hh  =  (tmax-t)/(dble(nstp*nskp)) 
nxt  =  0 

call  haming (nxt) 
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c 

if(nxt  .eq.  0)  stop  77 
write  (13,*)  ’haming  initialized’ 
c  integration  loop 
c 

do  1000  istp  =  l.nstp 
c 

c 

c  do  the  integration 

c 

do  500  jskp  =  l.nskp 
c 

call  haming(nxt) 
c 
c 

500  continue 

c 

c  PRODUCTION  CODE 
c 

c  normalize  variation 

dmag=dsqrt(x(3,nxt)*x(3,nxt)  +  x(4,nxt)*x(4,nxt) ) 
c 

do  600  i=l,2 

x(i+2,nxt)=  x(i+2,nxt)/dmag 
600  continue 

c 
c 

9  format(2x,5(e20.13,lx)) 

write  (1,9)  t,x(l,nxt) ,x(2,nxt) 
write  (2,9)  t,x(3,nxt) 
write  (3,9)  t,x(4,nxt) 
c 

c  move  state  vector  and  variation  vector  into  slot  1 

c 

do  900  i=l,4 
x(i , l)=x(i ,nxt) 

900  continue 
c 

c  reinitialize  haming 

nxt  =  0 

call  haming (nxt) 
c 

if (nxt  .eq.  0)  stop  77 
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write(13,*)  t,  ’haming  initialized’ 
c 

1000  continue 
c 
c 
c 

stop 

end 

include  ’haming_vp.for’ 
include  ’fftrhs_vp.for’ 
include  ’amat.vp.for’ 
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c 

c 

subroutine  rhs(nxt) 
c 

c  Van  der  Pol’s  equation  eom  and  equations  of  variation 
c 

common  /ham/t,x(10,4) ,f (10,4) ,err(10) ,nn,hh,mode,e 
double  precision  t,x,f ,err,hh,e 

common  /eval/  eig(2,4) ,npi(2) ,discl(4) ,dtbif ,iswap(2) 
complex* 16  eig.wp 

double  precision  disci ,dtbif ,wpp(2) 
common  /debug/  idebug 
double  precision  a(2,2) 
double  precision  xx(2) 
equivalence  (wp.wpp) 
c 

c  extract  state  and  variation 
c 

do  1  i  =  1,4 

xx(i)  =  x(i,nxt) 

1  continue 
c 

c  Van  der  Pohl’s  eom 

c 

f(l,nxt)  =  x(2,nxt) 

f (2,nxt)  ®-x(l ,nxt)+(e-e*x(l ,nxt)*x(l,nxt))*x(2,nxt) 
c 

c  A(t)  matrix 

c 

call  amat(nxt,a) 
c 

c  delta_x  dot  =  a  delta.x 

c 

c 

do  50  i  =  1,2 

f(i+2,nxt)  =  O.dO 
do  50  k  =  1,2 

f(i+2,nxt)  =  f(i+2,nxt)  +  a(i,k)*x(k+2,nxt) 

50  continue 
c 
c 
c 

return 

end 
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c 

c 

subroutine  amat(nxt,a) 
c 

c  calculate  A  matrix  for  Van  der  Pol’s  ’simple  ’  system, 
c 

common  /ham/  t,x(10,4) ,f (10,4) ,err(10) ,nn,hh,mode,e 
double  precision  t,x,f ,err,hh,e 
double  precision  xx(2),a(2,2) 
c 
c 

c  extract  state 
c 

do  5  i  =  1,2 

xx(i)  =  x(i,nxt) 

5  continue 
c 
c 

a(l,l)=0.d0 

a(l,2)=l.d0 

a(2,l)=-l.d0-2.d0*e*xx(l)*xx(2) 

a(2,2)=e-e*xx(l)*xx(l) 

c 

return 

end 
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c 


c 

program  pwr 
c 

c  do  power  spectral  density  calcs  on  deltax  data 
c 

dimension  data(20000) ,psd(10000) 
c 

c  read  in  data  file 

c 

read  (*,*)  npts.tmax 
n  =  npts/2 
do  100  i  =  l.npts 

read  (*,*)  t,data(i) 

100  continue 
c 

call  power (  data,  n,  psd,  tmax,  fmax  ) 
c 

c  output 

c 

open ( 1 , FILE= ’ psd . out ’ , STATUS= ’ UNKNOWN ’ ) 
df  =  f max/real (n) 
c 

do  200  i  =  l,n+l 

f  =  real(i-l)*df 
write  (1,1)  f,psd(i) 

1  format(lx,el2.5,  lx,  el2.5) 

200  continue 
c 

stop 

end 

c$INCLUDE:  ’power. for’ 
c$INCLUDE:  ’realft.for’ 
c$INCLUDE:  ’fourl.for’ 

include  ’power. for’ 
include  ’realft.for’ 
include  ’fourl.for’ 
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subroutine  powerCdata,  n,  psd,  tmax,  fmax) 
c 

c  one  sided  power  spectral  density  per  unit  time  of  real  function 

c  (2n  real  values)  stored  in  array  data.  Performs  fft 

c  on  real  function  using  realft  and  fourl,  then  returns 

c  power  spectral  density  in  psd,  n  values.  For  power  spectral 

c  density  PER  UNIT  TIME,  its  divided  by  tmax,  the  length 

c  of  the  time  interval.  Also  returns  fmax,  maximum  frequency, 

c  n  MUST  be  a  power  of  2!  Returns  n+1  values  of  psd,  starting 

c  at  frequency  0  and  going  to  fmax 

c  output  from  realft  is  in  cycles/sec,  converted  to 

c  RADIANS/SEC!,  psd  in  power  per  radian  per  unit  time! 

c 

dimension  data(2) ,  psd(2) 
c 

c  fft  of  real  function 
c 

call  realft (  data,  n,  +1) 
c 

c  conversion  to  per  radian  per  unit  time 
c 

div  =  l./(  2.*3.1415C*tmax  ) 
c 

c  extract  first  and  last  (real  valued)  transformed  points 
c 

psd(l)  =  data(l)*data(l)*div 
psd(n+l)  =  data(2)*data(2)*div 
c 

c  remaining  values  are  complex,  take  their  norm 
c 

do  100  i  =  2,n 
c 

i2ml  =  2*i-l 
i2  =  2*i 

psd(i)  =  (  data(i2ml)*data(i2ml)  +  data(i2)*data(i2)  ) 

1  *div 

c 

100  continue 
c 

dt  =  tmax  /  (2.*  real(n)  ) 
fmax  =  2. *3. 14159/  (2.*dt) 
c 

return 

end 


B-8 


c 

c 

subroutine  realft(  data,  n,  isign  ) 
c 

c  calculates  the  fourier  transform  of  a  set  of  2n  real  valued 
c  data  points.  Treuisliterated  from  c  version  of  numerical 
c  recipies. 

c  The  data  stored  in  dataCl . 2n)  are  replaced  by  the  positive 

c  frequency  half  of  the  complex  fourier  transform.  The  real 

c  valued  first  and  last  points  are  returned  in  data(l)  and 
c  data(2) .  n  MUST  be  a  power  of  2!  Also  calculates  the  inverse 

c  transform  of  a  complex  array  if  it  ’s  the  transform  of  real 

c  data,  but  the  result  must  be  divided  by  n. 
c 

c  isign  =  +1  for  time  ->  frequency  domain 

c  =  -1  for  frequency  ->  time  domain 

c 

dimension  data(2) 

c  double  precision  for  recurrences 


double  precision  wr,wi,wpr,wpi,wtemp, theta 
c 

cl  *  0.5 
c 

c  initialize  recurrences 
c 

theta  =  3.141592653589793d0/dble(n) 
c 

if (isign  .eq.  1)  then 
c  forward  transform 

c2  =  -0.5 

c  do  forward  transform 

call  fourKdata,  n,  +1) 

else 

c  set  up  for  an  inverse  ticuisform 

c2  =  0.5 
theta  =  -theta 

endif 

c 

c  prepare  to  separate  two  separate  transforms  of  alternating 

c  real  data,  auid  recombine  into  one  transform  on  positive  freq. 

c 

wtemp  =  dsin(  0.5d0*theta  ) 

wpr  =  -2.d0*wtemp*wtemp 

wpi  =  dsin(theta) 

wr  =  1 . dO  +  wpr 
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wi  =  wpi 
n2p3  =  2*n  +  3 
c 

c  deconvolution  loop,  case  i=l  done  separately  below 
c 

no2  =  n/2 
do  100  i  =  2,no2 

11  =  i  +  i  -  1 

12  =  1  +  il 

13  =  n2p3  -  i2 

14  =  1  +  i3 

c  the  two  separate  transforms  are  separated  out  of  data 

hlr  =  cl*(  data(il)  +  data(i3)  ) 
hli  =  cl*(  data(i2)  -  data(i4)  ) 
h2r  =  -c2*(  data(i2)  +  data(i4)  ) 
h2i  =  c2*(  data(il)  -  data(i3)  ) 

c  then  recombined  to  form  the  transform  of  the  original  data 

data(il)  =  hlr  +  wr*h2r  -  wi*h2i 

data(i2)  =  hli  +  wr*h2i  +  wi*h2r 

data(i3)  =  hlr  -  wr*h2r  +  wi*h2i 

data(i4)  =  -hli  +  wr*h2i  +  wi*h2r 

c  trig  recurrence 

wtemp  =  wr 

wr  =  wtemp*wpr  -  wi*wpi  +  wr 
wi  =  wi*wpr  +  wtemp*wpi  +  wi 
100  continue 
c 

c  final  fix  on  forward  transform,  or  inverse  transform 
c 

ifCisign  .eq.  1)  then 

c  squeeze  first  and  last  pts  into  first  data  slot 

hlr  =  data(l) 
data(l)  =  hlr  +  data(2) 
data(2)  =  hlr  -  data(2) 

else 

c  do  the  inverse  transform 

hlr  =  data(l) 

data(l)  =  cl+(  hlr  +  data(2)  ) 
data(2)  =  cl*(  hlr  -  data(2)  ) 
call  f our 1 (data,  n,  -1) 

endif 

c 

return 

end 
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c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


SUBROUTINE  FOURl(DATA,NN,ISIGN) 

fast  fourier  transform  for  complex  single  precision  data 
from  numerical  recipies. 

data:  input  real  vector  of  complex  valued  function, 
stored  in  time  ascending  order,  with  real/imag 
parts  alternating. 

nn;  the  number  of  complex  data  points,  data  will  actually 
have  2  nn  entries 

IT  IS  VERY  DESIRABLE  THAT  nn  BE  A  POWER  OF  2! 
isign:  direction  of  transform.  If  isign  =  +1,  goes  from 
time  ->  frequency  domain.  If  isign  =  -1,  from 
frequency  ->  time  domani  EXCEPT  output  transform  needs 
to  be  divided  by  nn  to  be  normalized. 

double  precision  for  trig  recurrences  only 
REAL*8  WR,WI,WPR,WPI,WTEMP, THETA 
DIMENSION  DATA(*) 

bit  reversal  section  of  routine 


N=2*NN 

J=1 

DO  11  1=1, N, 2 

IF(J.GT.I)THEN 

c  exchange  the  two  complex  numbers 

TEMPR=DATA(J) 

TEMPI=DATA(J+1) 

DATA(J)=DATA(I) 

DATA(J+1)=DATA(I+1) 

DATA(I)=TEMPR 

DATA(I+1)=TEHPI 

ENDIF 

M=N/2 

1  IF  ((M.GE.2) .AND.(J.GT.M))  THEN 

J=J-M 
M=M/2 
GO  TO  1 

ENDIF 

J=J+M 

11  CONTINUE 


c 

c  here  begins  the  Danielson-Lanczos  section  of  the  routine 
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recurrent  FFT’s  on  power  of  2  data  subsets 
MMAX=2 


outer  loop  executed  log  2  nn  times 


IF  (N.GT.MMAX)  THEN 
ISTEP=2*MMAX 

initialize  trig  recurrences 
THETA=6 . 283185307 17959D0/ (ISIGN*MMAX) 
WPR=-2 . DO+DSIN (0 . 5D0*THETA) **2 
WPI=DSIN (THETA) 

WR=1.D0 

WI=0.D0 


two  nested  inner  loops 

DO  13  M=1,HMAX,2 

DO  12  I=M,N,ISTEP 
J=I+MMAX 

the  Danielson-Lanczos  formula 
TEMPR=SNGL (WR) *DATA ( J) -SNGL (WI ) *DATA ( J+ 1 ) 
TEMPI=SNGL (WR) ♦DATA ( J+ 1 ) +SNGL (WI) ♦DATA ( J ) 
DATA ( J ) =D ATA ( I ) -TEMPR 
DATA ( J+ 1 ) =DATA ( 1+ 1 ) -TEMPI 
D ATA ( I ) =DATA ( I ) +TEMPR 
DATA ( 1+ 1 ) =DATA ( 1+ 1 ) +TEMPI 
CONTINUE 
trig  recurrences 
WTEMP=WR 

WR=WR^WPR-WI^WPI+WR 
WI=WI^WPR+WTEMP^WPI+WI 
CONTINUE 
MMAX=ISTEP 
end  outer  loop 
GO  TO  2 

ENDIF 


RETURN 

END 
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It  is  shown  that  the  concept  of  the  Lyapunov  exponent  can  be  extended  to  include  an  imaginary  part.  A  numerical  technique 
used  to  calculate  these  extended  Lyapunov  exponents  fm  second  order  systems  is  presented.  There  are  two  requirements  for 
this  extension:  the  definition  of  a  coordinate  frame  on  the  tangent  space  of  the  differential  equation,  and  an  extension  of  the 
cbssical  limit,  called  the  limit  of  the  mean.  An  application  of  the  technique  to  the  van  der  Pol  equation  for  the  constant 
coefficient  and  periodic  coefficient  cases  is  given.  The  extended  Lyapunov  exponents  found  using  this  technique  totally  agree 
with  the  eigenvalues  for  the  constant  coefficient  case.  In  the  peiiodic  coefficient  case,  not  only  do  the  extended  Lyapunov 
exponents  agree  with  the  Poincar6  exponents  calculated  using  standard  Floquet  theory,  they  confirm  that  the  imaginary  parts  of 
the  Poincar6  exponents  are  equal  to  the  quotient  ±i  where  r  is  the  period  of  the  trajectory.  This  imaginary  part  is  uncertain 
when  using  standard  Floquet  theory.  Additionally,  fast  Fourier  transform  (FFT)  techniques  are  used  to  validate  the  existence  of 
the  extended  Lyapunov  exponent  and  the  values  obtained  for  its  imaginary  part.  These  techniques  show  that  the  power  spectrum 
of  relative  motion  is  discrete  for  the  trial  cases  presented,  with  the  fundamental  frequency  almost  exactly  equal  to  the  calculated 
imaginary  part  of  the  extended  Lyapunov  exponent.  Coupled  with  the  successful  comparison  of  characteristic  exponents  for 
the  constant  coefficient  and  periodic  coefficient  cases,  this  power  spectrum  serves  to  decisively  validate  the  existence  of  the 
extended  Lyapunov  exponent. 
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